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The second-order solution of the problem of the inter- 
action of a train or regular waves with a completely sub- 
merged, horizontal circular cylinder in finite depth water 
is presented for the two-dimensional case. The incident 
wave is developed as a second-order Stokes 1 wave by use of 
a perturbation method and the solution of both the first- 
order and second-order scattering potentials is obtained 
numerically using the Green's function approach. The hydro- 
dynamic pressure and resulting first-order and second-order 
force coefficients are determined numerically and presented 
for various values of water depth, cylinder depth of 
submergence, and wave length. 
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SYMBOL INDEX 



Symbol Description 



a 


characteristic length, cylinder radius 


a 


dimensionless wavelength parameter, a = 2na/L 


b 


dimensionless constant, H = eb 


c i 


th 

dimensionless wave force coefficient in the i 
direction 


d 


dimensionless relative cylinder depth 


d 


cylinder depth 


i — 1 

u 


differential dimensionless arc length on the 
cylinder surface 


dC 2 


differential dimensionless length on the free 
surface 


f l 


first-order source strength function 


f 2 


second-order source strength function 


f* 


free surface particular solution portion of the 
second-order problem source strength function 


F( ) 


function of the quantities in parentheses 


F li 


dimensionless first-order wave force coefficient 
th 

in the i direction 


F 2i 


dimensionless second-order periodic wave force 

th 

coefficient in the i direction 


p ss 

F 2i 


dimensionless second-order steady-state wave 

th 

force coefficient in the i direction 


g 


acceleration of gravity 


G 


Green's function 


G* 


modified Green's function for the particular 
solution portion of the second-order problem 
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dimensionless relative mean water depth 
mean water depth 

first-order non-homogeneous boundary condition 

th 

function at the i nodal point 

dimensionless relative wave height, H = H/2a 

elevation difference between the wave crest 
and trough 

i, 

complex plane portion, i = (-1) 2 
wave number, k = 2 tt/L 

second-order non-homogeneous boundary condition 
function at the i*'* 1 nodal point 

dimensionless Bernoulli constant 

Bernoulli constant 

second-order dimensionless Bernoulli constant 
function of the quantities in parentheses 
dimensionless wave length 
wave length 

number of cylinder surface increments and 
nodal points 

number of free surface increments and nodal 
points 

unit normal vector on the cylinder surface in 
the outward direction 

spacial component unit normals in the horizontal 
and vertical directions, respectively 

order of 

dimensionless pressure coefficient 

Pressure 

principal value 



Symbol 



q 



r 

r • 



R 



e 



S 



1 



S 



2 



t 



t 



u 



1 



u 



2 



u 



2 



U 



x,y 



x,y 



a 

3 

A 

V 

6 

e 



S,n 






Description 
fluid velocity vector 
dimensionless linear distance 
dimensionless image linear distance 
real part of a complex number 
cylinder surface 
free surface 

dimensionless time, t = at 
time 

first-order complex potential function 

second-order periodic complex potential 
function 

second-order time independent complex 
potential function 

non-periodic function 

dimensionless special variables in the 
horizontal and vertical directions, respectively 

spacial variables in the horizontal and 
vertical directions, respectively 

complex matrix 

complex matrix 

increment 

vector operator 

phase shift angle 

perturbation parameter 

dimensionless spacial variables corresponding 
in direction to x and y, respectively 

dimensionless free surface elevation 
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Symbol 


Description 


n 


free surface elevation 


0 


plane polar angle 


X 


half surface interval 


y 


dummy of integration variable 




positive roots of y^.tan(y^h) - v = 0 


V 


dimensionless wave frequency 


7T 


constant, 3.14159 


p 


fluid density 


z 


summation of terms 


O 


wave frequency 


<P 


dimensionless velocity potential 


$ 


velocity potential 


Subscripts 

i 


direction of a force or pressure component 
or nodal point location, depending on context 


j 


nodal point location 


n 


normal ' partial derivative 


t 


time derivative 


x,x 


special partial derivative in the horizontal 
direction 


Y/Y 


spacial partial derivative in the vertical 
direction 


1 


first-order component or x spacial direction, 
depending on context 


2 


second-order component or y spacial direction, 
depending on context 
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Symbol 


Description 


Superscripts 

I 


incident wave component 


o 


homogeneous solution portion of second 
order scattering problem 


S 


scattering wave component 


* 


particular solution portion of second- 
order scattering problem 
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I. 



INTRODUCTION 



The solution to the two-dimensional linear wave/structure 
interaction problem has been well-studied for both the finite 
and infinite depth cases. The fluid motion resulting from 
the interaction of a train of regular waves with a submerged 
horizontal circular cylinder in infinite depth water v/as 
first studied by Dean [Ref. 1] . Ursell [Ref. 7] later studied 
the problem anew placing the solution on a rigorous basis. 

More recently, Ogilvie [Ref. 6] reconsidered the problem. He 
computed the first-order oscillatory forces and second-order 
steady-state forces for the following cases: (a) the cylinder 

held fixed, (b) the cylinder forced to oscillate sinusoidally 
in still water, and (c) the neutrally buoyant cylinder allowed 
to respond to wave motion. 

It can easily be shown that the steady-state part of the 
second-order forces can be computed from the first-order 
potential alone. Accordingly, Ogilvie did not obtain a 
complete second-order solution; the steady-state forces 
arise from the first-order velocity squared terms in Ber- 
noulli's equation. In addition to the steady-state forces, 
second-order oscillatory forces also exist which have a 
frequency twice that of the first-order forces and represent 
the subject of the present work. 

This thesis reconsiders again the submerged horizontal 
cylinder problem but the extension is made to include water 
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of finite depth. The cylinder is considered to be completely 
submerged and held fixed in a train of regular waves. The 
objective is to determine the complete second-order solution 
and, accordingly, determine the oscillatory second-order 
forces as well as the steady-state second-order forces. 

The problem is treated as a regular perturbation problem 
in the small parameter, incident wave height/ cylinder radius. 
Proceeding in this way the incident wave appears as a second- 
order Stoke 's wave. The boundary-value problems for both 
the first-order and second-order potentials are established 
and the solutions to both are obtained by use of the Green's 
function method. 
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II. THEORETICAL DEVELOPMENT 



A. FORMULATION OF THE PROBLEM 

A systematic representation of the problem considered is 
illustrated in Fig. (1) . A rigid right cylinder of radius a, 
submerged to a depth d in water of depth h, is subjected to 
a train of regular waves propagating in the positive x direc- 
tion. The basic problem is that of calculating the hydro- 
dynamic pressure on the immersed cylinder and the resulting 
force correct to the second-order in wave height of the 
incident wave. 

Assuming the fluid to be irrotational , a velocity 
potential, <J>, may be defined as: 

q = V<Hx,y,t) (1) 

/v 

where q denotes the fluid velocity vector. (The barred 
quantities denote dimensional quantities.) Moreover, assuming 
the fluid to be incompressible, the velocity potential must 
satisfy the Laplace equation, 

V 2 $(x,y,t) = 0 (2) 

within the fluid region. 

In addition to the differential equation, Eq. (2) , $ 
must satisfy certain boundary conditions. In specific, 
these are: 
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Figure 1: DEFINITION SKETCH 



1 . 



The z.ero normal velocity condition on the bottom, 
defined by y = - h, where FT denotes the mean fluid 
depth . 

2. The zero normal velocity condition on the surface, 
Sl (x,y) = 0, of the cylinder. 

3. The kinematic boundary condition on the free surface, 
S 2 (x,y,t) = n(x,t) = 0. 

4. The dynamic boundary condition on the free surface. 
These boundary conditions may be stated mathematically as 
follows: 



$-(x,-h,t) = 0 (3) 

V4>'VS 1 = 0 (4) 

9S 2 

q * v s 2 + = 0 (5) 

$^(x,Ti,t) + j[4>^(x,n,t) 2 + 4>-(x,n,t) 2 ] + gn(x,t) = gK 

( 6 ) 

where n denotes the elevation of the free surface, g denotes 
the acceleration of gravity, and K denotes the Bernoulli 
constant . 

For convenience in carrying out the solution, the 
variables are next made dimensionless using the cylinder 
radius, and the wave frequency, a, as follows: 
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x = x/a 


y = y/a 


d = d/a 


h = h/a 




H = H/2a 


K = K/a 


2-/ 

v = o a/g 


t = at 


(7) 


<J> = a$/ga 


0 = n/a 









H denotes the dimensionless wave height, K denotes the 
dimensionless Bernoulli constant, and t denotes the 
dimensionless time. 

Utilizing the parameters defined in Eq. (7), Eqs . (2-6) 

may be rewritten to concisely define the boundary-value 
problem in dimensionless form as: 



2 

V <Jt(x,y,t) = 0 in the fluid region (8) 

<j> (x,-h,t) = 0 (9) 

4> n (x,y,t) =0 on S 1 (x,y) = 0 (10) 

<J> x (x,n,t) n x (x,t) - <j> (x,n»t) + vn t (x,t) = o (li) 



<J> t (x, n / 1) + ^[<{> x (x,n ,t) 2 + <j> y (x,rw t) 2 ] + r)(x,t) = K (12) 

B. PERTURBATION EXPANSION 

Expanding <j> and ri in terms of a perturbation parameter, 
e, provides a means for converting the nonlinear boundary-value 
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problem into a series of linear problems. The solution to 
each linear problem may be tractable whereas the nonlinear 
problem is unsolvable. 

Since <p, ri , and K are functions of the small parameter, 
e, they may be written in the power series: 



<Mx,y,t;e) = J2 e n <{> (x,y , t) (13) 

n=l 



n(x,t;e) = X £% (x,t) (14) 

n=l n 



and 



K -- 





(15) 



In Eqs. (11) and (12) , <J> contains n, and, therefore, e 
implicitly. This can be converted to an explicit form in 
e by use of the Taylor series expansion: 



<J> (x,n,t) 



00 , , . m 

n (x,t ; g) 

m=0 ml 




4> (x,y ,t) . 
m J 



y=0 



(16) 



The perturbation parameter, e, is related to the v;ave height 
in an, as yet, unknown manner. 

Equations (13-16) as well as appropriate derivatives 
similar to Eqs. (13-16) may now be substituted into the 
boundary-value problem given in Eqs. (8-12) to obtain a 
series of linear boundary- value problems for <J)^, < f^, etc. 
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That is, upon equating coefficients of like powers of e, 
the following problems for the first two terms in the 
expansion for <{> can be precipitated: 

First-order (e) : 



V 2 4> 1 (x,y ,t) = 0 




(17) 


4> ly (x,-h, t) = 0 




(18) 


$ln (x,y , t) =0 on 


S- L (x,y) = 0 


(19) 


4>iy(x,o,t) - vn lt (x,t) 


= 0 


(20) 


n 1 (x,t) + 4> lt (x , o , t) = 


0 


(21) 


2 

-order (e ) : 


V 2 <f> 2 (x,y,t) = 0 




(22) 


4> 2y (x,-h,t) = 0 




(23) 


♦ 2n (x,y,t) =0 on 


S 1 (x,y) = 0 


(24) 



4> 2y (x,0,t) - vn 2t (x,t) = 

(25) 

-n 1 (x,t) 4> lyy (x,o,t) + 4> lx (x,o,t) n lx (x,t) 
n 2 (x,t) + 4> 2 1 (x / o , t) = -n 1 (x / t)<{» 1 (x,o,t) - 

n lt (x,t) <{> ly (x,o,t) - 2 + <(> ly (x,o, t) 2 ] + k 2 

(26) 
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There is a striking similarity between the first-order 
and second-order problems; only the right hand side of the 
free surface boundary conditions differ. In the first-order 
problem the free surface boundary conditions are homogeneous 
whereas the second-order free surface boundary conditions 
are non-homogeneous functions of the first-order velocity 
potential, first-order free surface elevation, and their 
derivatives . 

The first-order potential function may be represented 
by a function which is periodic and, therefore, the complex 
potential, u^(x,y), is defined as; 

4> 1 (x,y,t) = abR e [iu 1 (x,y)e _lt ] (27) 

where R e denotes the real part, a = 2ira/L, L denotes the 
wave length, and b is an unknown real constant. 

C. THE FIRST-ORDER BOUNDARY VALUE PROBLEM 

Since the first-order boundary-value problem is linear, 
the potential <f>^ may be expressed as a sum: 

*1 = $1 + (28) 

I S 

where <f>^ denotes the incident wave potential and <f>^ denotes 

the scattering potential due to the presence of the cylinder. 

In terms of the time independent complex potentials, using 

Eqs. (27-28): 
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( 29 ) 



1 , S 

U 1 = u i + u 



Remembering that the incident wave potential represents 
only the incoming wave, with all the effects of the cylinder 
represented by the scattering potential, then the incident 
potential must satisfy the first-order boundary-value problem 
when no rigid body is present. Therefore, satisfies 
Eqs. (17-18) and (20-21) , which represents simply the first- 
order boundary-value problem for a train of regular waves. 

The solution to this well-known problem in terms of the 
complex potential is: 



I _ _ 1 cosh [a (y+h) ] iax 
u i a cosh (ah) e 



(30) 



The relationship between v and a is defined by 



2 - 

v = — — = a tanh (ah) 



(31) 



which, in dimensional form, using a = ka, and k = 2 tt/L is: 



0 = gk tanh (kh) 



(32) 



Since the solution for u^ is known, the boundary-value 

s 

problem for u 1 may now be established. Substituting Eqs. 
(30), (29) and (27) into the first-order problem given by 

Eqs. (17-21), and eliminating n 1 between Eqs. (20) and (21) 
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results in a boundary-value problem .for the first-order 
scattering potential as follows: 



V 2 u^(x,y) = 0 



( 33 ) 



u^ y (x,-h) = 0 



(34) 



u 



S 

In 



(x,y) 



cosh 1 (ah) [ n y sinh t a (y +h )l + in x cosh (a (y+h) j e lax 
on S 1 (x,y) = 0 (35) 

u^y (x, 0) - vu^(x,0) = 0 (36) 



where n^ and n denote the x- and y-components , respectively, 

A A A 

of the unit normal vector, n = in + in , directed outward 

' x J y 

into the fluid. 



D. SECOND-ORDER BOUNDARY-VALUE PROBLEM 

In view of the linearity of the second-order problem, 
the same procedure as used in the first-order problem may 
be applied and leads to the representation of the second-order 
potential as the sum: 



<t> 2 = $2 + <t>2 (37) 

where (j)* denotes the second-order incident wave potential 

S 

and $2 denotes the scattering potential . 
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Considering the case where there is no body present, 

S S 

there is no scattering w T ave and therefore, <J>^ = $2 = 0 . 
Additionally, the boundary conditions on the immersed cylinder, 
Eqs . (19) and (24) are not applicable. When the substitu- 

tions of <{)^ for cj>^ and 4 * for $2 are made into Eqs. (22) , 

(23), (25) and (26), along with Eqs. (20), (21), (27) and 
(31), and upon eliminating n-^ and n 2 between Eqs. (25) and 
(26) , the resulting boundary-value problem for the second- 
order incident wave potential is: 

V 2 4>2 (x,y ,t) = 0 (38) 

<J>2y (>^/ “ h , t) = 0 (39) 

<f>2 y (x,0,t) + v4>2 tt (x,0,t) = - |b 2 (a 2 -v 2 )R e [ie i2(ax_t ) ] 

(40) 

The periodic solution to the incident wave boundary-value 
problem specified by Eqs. (38-40) is known and may be written 
as : 

,1 3 ,2 „ r . I, . -i2t n , . t . 

4>2 = jsb vR e [iu 2 (x,y)e ] (41) 

where the second-order complex potential, u* , is given by: 

u*(x,y) = - Jj- - c - ° sh [2a l y+h) 1 e l2ax (42) 

z sinh ‘ (ah) 
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Eqs. (41-42) are familiar in wave theory and exactly the 
same form as that given by Ref. 4 for second-order Stokes' 
waves. 

Having developed a solution for the second-order incident 



are substituted into Eqs. (22-26), and the previously deter- 
mined solutions for the incident wave potentials, Eqs. (30) 
and (42) , are utilized. After applying the relationships of 
Eqs. (20), (21), (27), (29) and (41) and eliminating n 2 

between Eqs. (25) and (26), the resulting boundary-value 
problem for the second-order scattering potential becomes: 



. 1 

wave potential, <j> 2 , the problem for second-order scattering 

IS IS 

potential may be determined. <J>^ = <J>i + ^ an< 3 <J> 2 = ^2 + ^2 



V 2 <f|(x,y,t) = 0 



(43) 



^2y = 0 



(44) 






4 • ,4 , ~ 

smh (ah) 




(45) 



on S (x,y) = 0 




1SS IS c I S 

— u, u. + u.u, - 6u. u, (46) 
v ly lyy 1 lyy ly ly 




+ U(x) 



on y = 0 
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where U(x) is a non-periodic function generated by the 

substitution of <j>^ into Eqs . (25) and (26) . For brevity, 

this function is not written out since it is not needed in 

determining the second-order pressures and forces. 

g 

The boundary-value problem in given in Eqs. (43-46) 

-i2t 

is time dependent according to e , except for the U(x) 

g 

term in Eq. (46). Therefore, the solution for <j> 2 may be 
taken in the form: 

.S 3 . 2 _ . r S, , -i2t . ~S , .. 

<J> 2 = ^ ab vR 0 i [u 2 (x,y) e + u 2 (x,y)] (47) 

where the last term denotes the time independent portion of 

the complex potential. Separate boundary-value problems for 

s 

u 2 and u 2 arise from the substitution of Eq. (47) into Eqs. 
(43-46) . However, as will be demonstrated in the pressure 

g 

and force development, only the <j> 2 term is required. There- 

g 

fore, the time independent part of <}> 2 will be zero, and as such 

will not be further developed; concentration will be directed 

g 

towards the solution for u 2 . 

When Eq. (47) is substituted into Eqs. (43-46), the 
resulting boundary-value problem for the second-order 
scattering potential is: 

V 2 u 2 (x,y) = 0 (48) 

u 2y (x,“h) = 0 (49) 
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u? (x,y) = fn sinh [2a (y+h) ] 

2n sinh 4 (ah) L y 

+ in x cosh [2a (y+h) ]J e l2ax on S 1 (x,y)=0 (50) 

u 2y ( x ' 0) " 4vu^ (x , 0) = f*(x) (51) 



where 



, 2a r l S S IS . I S IS I .... 

f * ( x ) = — I — u -.U-, + u,u, - 6u, u, + — u , u (52) 

KJ 3v v ly lyy 1 lyy ly ly v ly lyy v ' 



, I S o S 
“ 4u lx u lx " 2u lx 



- 3u 



iy 



] y=0 



Again, there is similarity in form between the first-order 
and second-order problems; the only difference in form being 
that the second-order free surface boundary condition, 

Eq. (51), is non-homogeneous . 

Since Eq. (51) is non-homogeneous, further use of the 

S 

linear superposition theory is made by defining u 2 as the 
sum; 



S S° , S* .... 

u 2 = u 2 + u 2 (53) 

g® g * 

where u 2 denotes the homogeneous solution and u‘ 2 the 

particular solution to the boundary-value problem as stated 

S * s° 

in Eqs. (48-51). More precisely, u 2 and u 2 are defined as 
the solutions to the boundary-value problems obtained from 
the substitution of Eq. (53) into Eqs. (48-51) as follows; 
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s * 

Particular Solution (u^ ) : 



O Q * 

v U 2 (x ,y) = 0 


(54) 


U 2 * (x,-h) = 0 


(55) 


U® (x,0) - 4vu^| (x , 0 ) = f*(x) 


(56) 


Homogeneous Solution: 




2 S° 

V u* (x,y) = 0 


(57) 


s® 

u 2y (x,_h) = 0 


(58) 


s° s° 

U- 2 y(x, 0 ) - 4vUp (x,0) = 0 


(59) 



u^ (x,y) = i \n sinh [2a (y+h) ] 

sinh (ah) L y 

•k 

+ in^ cosh [2a (y+h) ] Je i2ax - ^ 2 n (x,y) (60) 

on S^(x,y) = 0 

By the judicious division into homogeneous and particular 

• S* 

solutions, the non-homogeneous problem for U 2 contains no 

boundary condition on the cylinder surface. This boundary- 

s * 

value problem for as stated in Eqs. (54-56) is identical 

to that associated with the linear problem for the potential 
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resulting from a harmonic pressure variation of amplitude 
distribution f* (x) on the free surface in water of depth h. 

The homogeneous boundary-value problem for u^ , defined by 
Eqs. (57-60) , is similar in form to the first-order problem 
given by Eqs. (33-36) , the significant difference being the 
term 4v in place of v which occurs in the free surface boun- 

g° 

dary condition. Thus, the method of solution for U£ will 
be similar to that of the first-order scattering potential. 

E. DEFINITION OF THE INCIDENT WAVE 

Having developed the first-order and second-order boun- 
dary-value problems, it is now appropriate to completely 
define the incident wave height and in the process determine 
expressions for the unknown constants b and • Solving 
Eqs. (21) and (26) for and respectively, and then 

substituting the results into the expression for the free 
surface elevation as given by Eq. (14) yields: 

n(x,t) = - e* lt + e 2 [-* 2t + $ lt * lyt + *itt*ly (61) 

- + + ly ) + K 2 ] + 0 (e 3 ) 

where the potentials are evaluated at y = 0 . Eq. (61) is an 
expression of the free surface elevation in terms of the 
total potential, and therefore, includes the effects of both 
the incident and the scattering waves. Hence, as it is of 
interest to obtain an expression for the free surface elevation 
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of the incident wave, the scattering potentials are set to 
S S 

zero, i.e. <J>^ = 4> 2 = 0 . Evaluating Eq. (61) using the known 
solutions for the incident wave potentials, Eqs . (30) and 

(42), along with Eqs. (27) and (41) , then yields: 



1 , .. r i(ax-t), , 2fb^(v^-a^) 

n (x,t) = ebRje ] + e [ — — - + 



K, 



(62) 



ab^ cosh (ah) [2 + cosh (2ah) ] D , i2(ax-t).‘l 
+ R [e J 

* sinh- 3 (ah) e J 



+ 0(e J ) 



where ^^(x^) denotes the dimensionless sxirface elevation 
for the incident wave with no cylinder present. If the x-axis 
is placed at the mean water line then the constant term must 
vanish and, therefore, 



K 2 = 



. 2 , 2 2 . 

b (a - v ) 

4 v 



(63) 



The remaining terms in Eq. (62) are time dependent periodic 
functions, the second-order term at twice the frequency of 
the first-order. 

Defining the dimensionless wave height given in Eq. (7) 
as the difference in elevation between the crest and trough 
of the incident wave, then eb must represent the wave ampli- 
tude. The second-order term in Eq. (62) at twice the funda- 
mental frequency makes equal contributions to surface elevation 
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at both the crest and trough, so that it contributes nothing 
to the wave height. Thus, in terms of dimensionless wave 
height, b is given by 



H = eb = H/2a 



(64) 



where H denotes the elevation difference between the wave 
crest and trough. 

Expressing the incident wave profile in terms of the 
dimensionless wave height yields: 



It may be noted that Eq. (65) agrees with the expression 
for the second-order Stokes' wave given by, for example, 

Ref. 4 . 

F. PRESSURES AND FORCES 

The pressure may be formulated by use of Bernoulli's 
equation, arranged as follows: 

1 2 2 _ _ 

P(x,y,t) = -P$— - =-p[$— + $— ] - pgy + pgK (66) 

x y 

where P(x,y,t) denotes the dimensional pressure and p denotes 
the fluid density. By use of Eqs. (7) , (13-16) , (27) , (41) , 



ri*(x,t) = H cos (ax-t) 



(65) 




[2 + cosh (ah) ] cos [2 (ax-t)] 
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(47), and (64), Eq. (66) may be reduced to dimensionless form, 
and carried to the second-order in wave height as follows: 



P(x,y,t) = -y - HaRju^ 



-it 



(67) 



4v 





In Eq. (67) p(x,y,t) denotes the dimensionless pressure 



coefficient and is defined as: 



p(x,y,t) 



_ P (x,y ,t) 
pga 



( 68 ) 



The first term in Eq. (67) represents the hydrostatic 
pressure as y is the dimensionless depth beneath the mean 
free surface (y = 0) . The second and third terms are harmonic, 
the third term having twice the frequency of the second and 
representing the harmonic second-order contribution. The 
remaining terms in Eq. (67) are independent of time and 
provide the time-average or steady-state force components. 

Expressing the components of the wave force vectors in 
terms of integrals of the pressure over the cylinder surface 
area, the dimensionless force coefficients may be written 
as: 




i = 1,2 



(69) 
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where the dimensionless force coefficients in the x and y 
directions are defined, respectively, as: 



C-^t) 



Vt> 

Pga 3 



C 2 (t) 



yt) 

pga 3 



(70) 



(71) 



with F x (t) and Fy(t) denoting the force components. Addi- 
tionally, dC^ denotes a dimensionless differential arc 
length along the circumference of the cylinder. 

Applying Eq. (67) to Eq. (69) yields: 



C ± (t) 





)e- i2t 



+ 






1 ] n . dC n + 0 (H 3 ) 

l 1 



(72) 



As the first term in Eq. (72) results from the hydrostatic 
pressure on the cylinder, it represents simply the buoyancy 
force. Omitting the hydrostatic force, the hydrodynamic 
force may be written as: 
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(73) 



C ± ( t) = H F li cos(6 li -t) + H 2 

+ H 2 [F 2i cos ( 62i-2t) + F®?] + 0 (H 3 ) 



where the first-order, second-order, and steady-state force 
coefficients and phase shift angles are defined by comparison 
of Eqs. (72) and (73) as follows: 



i 6. . r 

,.e 1 = a \n, n . dC n 

li ~ J 111 



(74) 



l6„ . 2 r r 2 

2i a /, 6v - \ 

; ± e = " Ul -- “ u-,„)n,dC 



6v" 2 2, 

— u 2 ' u lx ' u ly'“i““l 



(75) 



2 2 
V 



F 2i ‘ ft C / ( I“1 X I + Kyi + K “ 1)n i dc l 



(76) 



The dimensionless force coefficients, F^ and F 2 ^ are real. 

The first term in Eq. (73) represents the first-order 
forces of the fundamental frequency, a. The periodic portion 
of the second term in Eq. (73) represents the second-order 
contribution to the force at twice the fundamental frequency. 
The last second-order term represents the steady-state or 
time- independent contribution, sometimes called drift force. 

Evaluation of the force coefficients defined by Eqs. 
(74-76) is the main objective of this thesis. Therefore, 
it is necessary to evaluate the first-order and second-order 
complex potentials, u 1 and u 2 , as well as derivatives of u^ 
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on the surface of the cylinder. Additionally, evaluation 
of u^ and its derivatives on the mean free surface (y = 0) 
is required. 
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III. METHOD OF SOLUTION AND NUMERICAL PROCEDURES 



A. SOLUTION OF THE FIRST-ORDER PROBLEM 

The use of a Green's function to express the solution 
to the first-order boundary-value problem was first formu- 
lated by John [Ref. 5] , and applied to submerged ellipsoids 
by Garrison and Rao [Ref. 3] . This method is considered to 
be applicable, in principle, to arbitrary shapes and is 
mathematically the most straightforward . The Green's func- 
tion, G, of unit strength which satisfies Eq. (33) as well 
as the boundary conditions, Eqs. (34) and (36) is given by: 



00 r- 



Gts.yfS.niV) - In r - In r' + 2PV f Ieo shjji(gjin 

q*' Licosnyh) tvcoshyn-ysxnn yh; 



-yh sinh (pn) sinh(yy) 
sinlTIyhT 



H | 

- cos | x- 



£ I dy 



(77) 



- l 



4ir 



2a^h+sinh (2a^h) 



cosh [a^ (y+h) ] cosh [a.^ (n+h) ] cos [a.^ | x-£ | ] 



where : 



r 



r’ 



1 



[(x-O 2 + 


(y-n) 2 ] 2 


(78) 




l 




[(x-£) 2 + 


(y+n) 2 ] 2 


(79) 



The symbol, a^, is defined in terms of h and v as the 
solution to the equation: 
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F ( a 1 , h , v ) 



0 



(80) 



where F is defined by: 

F(a^,h,v) = tanh(a^h) - v (81) 

Comparison of Eqs. (31), (80) and (81) demonstrates that a^ 

is clearly the equivalent of a. However, this will not be 
the case for the second-order problem, requiring the use of 
separate notation. In Eq . (77), PV denotes the principal 

value of the integral. 

This Green's function was also given in series form by 
John [Fef. 5] as: 

« (y 2 +v 2 )cos[y (y+h) ] -y |x-£| 

G(x,y?£,n;v) = 2 tt 2^ * ^ cos[y (n+h)]e 

k=l vy k ~hy k ('^+v ) 

ia |x-£| 

4-rrcosh[a 1 (y+h) ]cosh[a, (n+h) ]e 

- i ± (82) 

2a^h + sinh(2a^h) 

where a^ is as defined in Eq. (80) and the quantities y^ 
are defined as the real positive roots of the equation: 

K(y k ,h,v) = 0 (83) 

where the function K is further defined as: 
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K(y,h,v) = y tan(yh) + v 



(84) 



g 

Following the Green's function method of solution, u^ 
is written as the integral over the cylinder arc length, 
C^, as: 



u^x^) = 



_ 1 _ 
2 TT 



A 



(£,n)G(x,y;£,n;v) dC 1 



(85) 



where (£,y) denotes points on the immersed surface, f^(£,n) 
denotes the unknown source strength, and dC^ = dC^/a denotes 
the differential arc length on the cylinder surface, made 
dimensionless with the characteristic length, a. 

Although the solution to the first-order boundary- value 
problem as stated in Eqs . (33-36) is given by Eq. (85) , the 

source strength function, f^(£,n), must be determined in 
order to evaluate the potential. From potential theory, 

S 

the normal derivative of the potential, u^(x,y) , on the 
surface of the cylinder is: 



u ln (x,y) = \ f A X,y) + Jif J f l (?' r >)G n (x,y;C,n;v)dC 1 (86) 

c ** 

where G , the normal derivative of G, may be determined by 
differentiation of either Eq. (77) or (82) in a straight- 
forward manner. 

Applying the final boundary condition, Eq. (35) , leads 
to an integral equation which may be solved for f^, 
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c 



ff L u,n) 



G n (x,y ; £, n ,v) d^ 



(87) 



1 



cosh(a ~ h) , [ n y sinh[a(y+h) ] + in x cosh ( a (Y +h ) lje iax 



The solution for from Eq . (87) may then be used in Eq . 



B. SOLUTION OF THE SECOND-ORDER PROBLEM 

The similarity in form of the boundary-value problem 
for the homogeneous part of the second-order potential, 

Eqs. (57-60), to the first-order potential is now utilized. 
Since the only significant difference occurs in Eq. (59) , 

S° 

where v is replaced by 4v, we may represent in a form 

similar to Eq. (85) as: 



where G (x,y ; £,n ; 4v) is defined by replacing v by 4v in Eqs. 
(77) and (82). Therefore, a 2 is defined by: 



and a^ is replaced by a 2 in Eqs. (77) and (82) also. In the 
case of Eq. (82) , is defined as the real positive roots of: 



S 

(85) to determine the potential, u^. 




1 



( 88 ) 



F(a 2 ,h,4v) = 0 



(89) 
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K(y k ,h,4v) 



0 



(90) 



The source strength function, f 2 (£,n), appearing in 
Eq. (88) is again determined by applying the kinematic 
boundary condition on the cylinder surface as given by Eq. 

(59) . The integral equation for f 2 may then be given by: 

j f 2 (x,y) + — J± 2 (£,n) G n (x,y; £,n;4v) dC^ (91) 

C 1 

"n^sinh [2a (y+h) ] + in^cosh [2a (y+h) ] 
smh (ah) 

on S (x,y) = 0 

S* 

However, to solve Eq. (91) for f 2 , first a solution for u 2 

s * 

must be obtained and u 2n evaluated on the cylinder surface. 

S* 

Thus, at this point the solution to u 2 is sought. 

s* 

The complex potential u 2 is defined as the solution to 
the boundary-value problem, Eqs . (54-56). The form of this 

problem is recognized as being equivalent to that of fluid 
motion produced by a periodic pressure distribution (as 
specified by f*(£)) with frequency 2a on the free surface. No 
cylinder is considered to exist in the fluid region. 

The solution to this problem may be formulated as an 
integral over the free surface extending from negative 
infinity to positive infinity. Using the present 



i2ax 



S * 

U 2n (x,y) 
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s* 

dimensionless representation, u 2 becomes: 

u 2* (x ' y) = ^*(5)G*(x,y;?/0;4v) d£ 



( 92 ) 



The Green's function, G* , for the problem is given in Ref. 8. 
Upon transforming the solution to the non-dimensional form: 



G*(x,y,5,0;4v) = -2PV / cosh [myth) ] cos [ u (x-q ]du 

-«■ 4vcosh(yh) - ysinh(yh) 



4ircosh [a 2 h] cosh [a 2 (y+h) ] cos [a 2 (x-£) ] 



-l 



2a 2 h + sinh(2a 2 h) 



(93) 



where a 2 is defined by: 



F(a 2 ,h,4v) = 0 



(94) 



The normal derivative of u 2 is required in order to 
completely define the kinematic boundary condition on the 
cylinder as given in Eq. (60). From Eq . (92), this becomes: 



U 2n (x,y) = Tii >(OC n( x, y ;? ,°,.4v) d£ 



(95) 



The derivative of G* in the direction normal to the cylinder 
surface may be obtained by differentiation of Eq. (93) in a 
straightforward manner. Thus using the values for f*(£), 
as defined in terms of the first-order solution only, as 
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s* s* 

given in Eqs. (52) and (56), and U 2 n may be determined, 

S * 

Using / U 2 may then be determined to complete the 

second-order boundary-value problem solution. 



C. NUMERICAL METHODS 

Determination of the forces on the cylinder surface 
requires solving for the potentials u^, U 2 and the derivatives 
of u-^ at points on the surface, as shown in Eq. (72). Thus, 
the problem is now one of solving for both the first-order 
and the second-order scattering potentials, which requires 
solving for the source strength functions f^ and as 

as the free surface pressure distribution source strength 
function, f*. In view of the complexity of the equations 

it is natural to attempt a numerical solution. 

S 

The first-order scattering potential, u-^ , as given by 
Eq. (85) is dependent upon the first-order source strength 
function, f^. Thus, it is necessary to solve Eq. (87) for 

S 

f-^ to determine u^. A numerical solution may be developed 
by dividing the cylinder surface into elements of length 
A6 = 2iT/m, with the center of each element assigned an index. 
Since f^(£,ri) is recognized as a well-behaved function for a 
smooth surface cylinder, it is appropriate to define the 
following: 



a ij (v) 



IT 



AC 



g (x. ,y . ; £,n;v) dc n 
J n 1 1 I 



(96) 



19 



i , j 1,2,3,. ..,m 
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(97) 



h i ° cosh ' (ah) ' [ n y (x i' y i ) sinh[a(y 1+ h)] 

iax. 

+ cosh [a (y^+h) ] e 1 

£ li " £ l (x i' y i> 

and thus Eq. (87) may be approximated by the complex matrix 
equation: 

oi.^j(v)f^j — ^^lj i / D — 1 , 2 , . . . ,m (99) 

Upon evaluation of ou^ (v) using Eq . (96), the inversion of 

Eq. (99) may be carried out on the digital computer to 
determine f^ at each nodal point on the surface of the 
cylinder . 

For the purpose of evaluating a and 6, either Eq. (77) 
or Eq . (82) may be used. For evaluation of a and 8 when i = j 

Eq. (77) must be used to allow separate treatment of the 
logarithmic singularity as r approaches 0. The difficulty 
encountered with the logarithmic singularity may be overcome 
by carrying out its integration analytically. Garrison 
[Ref. 2] showed for the calculation of that the contribu- 

tion of the normal derivative of the In (r) in Eq. (77) 
integrated over the singular element of AO is A6/2, not only 
for the nodal point (x^,y^) = (£,n)/ but for all nodal points 
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( x i,y±) on* the cylinder surface. Additionally, to calculate 



integral of the In (r) over the singular element of A0. 

A second singularity arises in the evaluation of the 
infinite integral occurring in Eq. (77) . The first term in 
the integrand is singular at p = a. Since the numerator 
approaches zero like (y-a) near a. Garrison [Ref. 2] sub- 
tracted l/(y-a) from the singular term in the range 0 <_ y <_ 2a 
and added the contribution of the integral of l/(y-a) which 
in this case was zero. This converted the singular integrand 
to a regular function. Applying this technique to numeri- 
cally evaluate the resulting integral between 0 and 2a, 
Simpson's three-eights rule is used with an odd number of 
subdivisions, thus ensuring that the point y = a is not 
encountered. 

With f^ determined from the inversion of Eq. (99) , the 
first-order potential and its derivatives may be evaluated 
on the cylinder surface. Replacing the surface integrals 
with summations, Eq. (85) and its derivatives may be written 
as : 





i , j 1/2,3 



f • • • / m 



( 100 ) 




( 101 ) 




( 102 ) 
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s s 

where u^, u^ etc. denote functions evaluated at the 
th 

i nodal point on the cylinder surface. The complex 
matrices in Eqs . (100-102) are defined as follows: 



hi™ 



e xij (v) -57 



AC lj 

y^x( x i ,Y i ; ?'0;v)dC 1 



if j if 2 



ac i 5 



B yij ™ - 57 AC , 

AL-ij 



^ r V) dC^ 



,m (103) 



(104) 



(105) 



The first-order incident potential and its derivatives may 
be evaluated directly at each nodal point after differentia- 
tion of Eq. (30) as needed, and combined with the scattering 
potential and its respective derivatives. 

To solve the second-order problem, the function f*(£) 
required in Eq. (95) and as defined by Eq. (52) must be 
evaluated numerically. Dividing the mean free surface 
(y = 0) in the vicinity of the cylinder into n equal incre- 
ments, f* is evaluated at each nodal point. The numerical 
solution uses again Eqs. (100-105) with y^ = 0 for the 

g 

purpose of evaluating u^ and its derivatives at the mean 
free surface nodal points. The incident first-order potential 
and its derivatives are evaluated directly at each nodal 
point, and combined with the scattering potential and its 
derivatives to solve Eq. (52) for f*. 
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Having determined f* at all nodal points on the free 

s * 

surface, the boundary-value problem for ^ given by Eqs. 

S* S* 

(54-56) may be solved, i.e. ^ and U 2 n may be evaluated 
on the cylinder by solution of Eqs. (92) and (95) respec- 
tively. Expressing the integrals as numerical summations in 
complex matrix form yields: 



u 2ni = j 



x 1,2,..., m 
j ■— 1,2,. ,.,n 



(106) 




f*e* . ( 4v) 
D ID 



where 



(107) 



* 

o i j (4 v) 
( 4 v ) 




Y i : C/n;4v) 



/Y i ; 4v) 









(108) 



(109) 



s* s* 

Thus, using f * , both U 2 and U 2 n may be directly evaluated 

at each cylinder surface nodal point. 

To solve the boundary-value problem for the second part 

S° 

of the second-order scattering potential, , as specified 

by Eqs. (57-60) , Eq. (88) must be evaluated. However, to 

S° 

evaluate Eq . (88) for ^ , the second-order source strength 

function, f 2 , must be evaluated by numerically solving 
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the integral equation, Eq . (91) . Replacing the integral 

equation with a complex matrix summation: 

f 2i + f 2j“ij^ 4v ^ = 2k i i/3 = (110) 

where cuj (4v) is defined by Eq. (96) , 



f 2j = f 2 < x j'Yj> 



(HD 



k. = 
i 



j ("n (x. ,y . ) sinh [2a(y .+h) ] 

4 ^hi Li 11 1 



( 112 ) 



sinh (ah) 



i2ax . 



1 g * 

+ in x (x i ,y i )cosh[2a(y i +h) ]J e - u 2n (x i ,y i ) 



Once a^j(4v) is evaluated, Eq. (110) may be inverted to obtain 
the second-order source strength function, f ^ < at each nodal 
point on the cylinder surface. Stating Eq. (88) in summation 
form: 




f 2j 6ij(4v) 



i , j 1,2,3, ...,ro 



(113) 



where 6. . (4v) is defined by Eq. (103) . Using Eq. (113) , 

S® 

U£ may be determined at each nodal point on the cylinder 

s * 

surface, then combined with u^ and the second-order inci- 
dent potential, u^ , to evaluate the total second-order 
potential, ^ • 
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Since the first-order potential, u^ , and its derivatives, 

and u ^y/ as well as the second-order potential, are 

now known at each nodal point on the surface, the first- 
order, .second-order and steady-state force coefficients and 
phase shift angles may be determined using Eqs. (74-76) . As 
shown earlier, the integrals may be replaced by summations, 
using nodal point values and the arc length increment, A0. 

Once each integral is evaluated, then the force coefficient 
becomes the absolute value or modulus of the complex integral 
result and the phase shift angle is the angle whose tangent 
is the imaginary part over the real part of the complex 
integral result. 

D. COMPUTER SOLUTION 

In order to utilize the method of solution described 
in Sections B and C, a computer program was developed to 
carry out the indicated calculations. Although accuracy was 
a prime consideration, an equally important requirement was 
to minimize computer run time. To achieve this, every attempt 
was made to utilize symmetry and to generate certain constants 
and matrices to eliminate redundant computations. 

The program consists of a main program that flows in the 
order of computation presented in Sections B and C, along 
with four subroutines. Two subroutines to solve the first- 
order Green's function, using both forms of the function, 

Eqs. (77) and (82), were developed by Garrison [Ref. 2], 

These subroutines, GREEN and GREENS respectively, have been 
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modified to include calculation of the first-order scattering 
potential derivatives, evaluation of the first-order scattering 
potential and its derivatives on the mean free surface, y = 0, 
and evaluation of the modified Green's function, G*, for 

S* g ★ 

the determination of u 2 and u 2n at each cylinder nodal 
point. 

For elements of the a and 3 matrices corresponding to 
small values of the parameter (x~£) , GREEN is used while for 
larger values of (x-£) , GREENS is used. With the exception 
of the diagonal elements in Eqs. (96) and (103-105) and a few 
surface nodal points in Eqs. (108-109), the majority of 
elements of a and 8 are calculated by GREENS. This is most 
fortunate as the series form converges rapidly, requiring much 
less computer time than the integral form. 

Subroutine GEODAT reads the input geometrical data, 
generates the matrices hu and as defined by Eqs. (97) and 
(112) respectively, and calculates certain geometrical 
parameters and matrices for repeated use in GREEN and GREENS. 
Subroutine COMAT inverts the complex matrix equations and 
thus is used to invert Eqs. (99) and (110) to determine 
f^ and f 2 respectively. 

A cross-reference between text and computer program 
nomenclature is given in Table I. 
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TABLE I: Computer Program — Text Symbol Cross-Reference 



Text 


Computer 

Program 


Text 


Computer 

Program 


a 


A 


n 


ANX(I) 






X 


d 


D 


n 


ANY (I) 






y 


f l 


F(I,1) 


U 1 


U1(I) 


f 2 


Fl (1,1) 


U lx 


U1X (I) 


f* 


FS (L) 


u ly 


U1Y (I) 


F 11 


Cl(l) 


u* (surface) 


U1IS 


F 12 


Cl (2 ) 


u lx ( sur ^ ace ) 


U1ISX 


F 21 


C2 (1) 


I 

Ufy (surface) 


UlISY 


F 22 


C2 (2) 


u 3yy ( Surface - 


U1ISYY 


F SS 

21 


C3 ( 1) 


g 

u^ (surface) 


U1SS 


SS 

F 22 


C3 (2) 


s 

Ufx (surface) 


U1SSX 


G 


GI J , GI JEXT 


s 

u ly ( Sur ^ ace ) 


U1SSY 


G* 


GIJ , GI JEXT 


s 

u iyy ^ surface ^ 


U1SSYY 


G n 


GNIJ , GNJI 


u 2 


U2 (I) 


G 

X 


GXI J , GXJI 


X 


X(I) 


G y 


GYI J , GY JI 


y 


Y (I) 


G yy 


GYY 


a 


ALPHA ( I, J) 


h 


H 


3 


BETA (I , J) 


2hi 


HH (1,1) 




BETAX ( I , J) 


k i 


PK (1,1) 


e y 


BETAY (I, J) 


m 


NPTS 


A0 (cylinder) 


DELTHE 


n 


NSPTS 


hE, (surface) 


DELX 
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Text 


Computer 

Program 


Text 


Computer 

Program 
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IV. DISCUSSION AND RESULTS 



A. SELECTION OF COMPUTER PROGRAM PARAMETERS 

The computer program was developed on the IBM- 360 
computer using FORTRAN H. All subroutines were compiled on 
DATA CELL to minimize compile time and core size. In order 
to maximize accuracy while minimizing computational time, 
two key computer parameters were evaluated to determine 
optimum values; cylinder nodal points and free surface model 
points . 

Using the first-order solution only, the number of nodal 
points on the cylinder surface was varied from 12 to 60 and 
compared with the results determined by Garrison [Ref. 2]. 
The minimum number- of cylinder surface nodal points that 
provided accuracy within one percent of- those obtained by 
Garrison using 60 points, was 24. Accordingly, the good 
correlation from 60 points down to 24 points led to the 
selection of 24 nodal points on the cylinder surface. 

Since the free surface integral was to be carried out 
from -oo to +«>, the next step was to establish the finite 
interval of convergence on the surface as well as the sub- 
division size. After the outer limits of the free surface 
integral were determined, the subdivision size was varied 
from 4 to 128 subdivisions per second-order wave length 
holding the total interval constant. Above 64 subdivisions 
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the results did not vary more than one percent, leading to 
the selection of 64 subdivisions per second-order wave length 
oh the free surface. 

Additionally, the total interval, 2X, was varied from 
one to eight wave lengths in both the positive and negative x 
directions. The second-order force coefficients and their 
respective phase shift angles converged to a small value 
varying periodically. Convergence to this limit cycle 
occurred between the first and second wave length out from 
the center of the cylinder. Therefore, a total interval of 
from -2 to +2 second-order wave lengths was chosen as adequate 
for generation of the numerical results. 

Although the values of the second-order force coefficients 
and their respective phase shift angles tended to converge 
with increasing the limits of the free surface integration, 
there remained, in general, the small limit cycle as noted 
above. The amplitude of the cycle was a function of the 
parameter, a, and was of significant magnitude only for 
values of a less than 0.25. The effects of a on the limit 
cycle for one of the second-order parameters, horizontal force 
coefficient, is demonstrated in Fig. (2) . This represents a 
plot of percent variation of the force coefficient from the 
mean versus the surface interval size, X/(L/2), corresponding 
to values of a of 0.25, 0.20 and 0.15. Figure (2) is 
representative of the behavior of all second-order force 
coefficients and phase shift angles and demonstrates the 



54 



0.15 



a x 





vo 



— If\ 



3 

.5 



— fA 



- OJ 



4 

> 



EH 

£5 



W 

O 

s 

g 

CO 



LT\ 



lf\ 

I 



O 

r— • 

I 



Lf\ 

T— 

I 



NV3N W0H3 NOIiVIAaa 1113033 d 



55 



Figure 2: SECOND-ORDER HORIZONTAL WAVE FORCE COEFFICIENT VERSUS SURFACE INTERVAL 



increase in the amplitude of the variation with a decrease 
in a. The amplitude of the periodic variation ranges from 
up to fifty percent at a wave number of 0.15 and depth of 
d = 1.4, to less than two percent for all values of wave 
number above 0.25. 

B. RANGE OF APPLICABILITY 

In a complex computer program of the type developed to 
carry out the present numerical scheme, certain limits with 
respect to numerical stability, etc. show up. These computing 
limits arise for various reasons such as overflows caused by 
computing hyperbolic functions of very large numbers, 
numerical convergence instabilities, etc. While no attempt 

was made to investigate each of these numerical problems, a 

\ 

list of the limitations of the particular program developed 
in this thesis to carry out the computations are as follows: 

minimum a — dependent upon acceptable error, but 0.25 
or less 

maximum a — deep water condition reached 
minimum h — 3 

maximum h — deep water condition reached for a > 0.25, 
h = 20 for a < 0.25 

minimum d — from 1.4 at h equal to or less than 5 down 
to 2 . 5 at h = 20 

minimum SMIN — 0.12 when cylinder depth is other than 

near the free surface 

maximum SMIN — 0.30 when cylinder depth is near the 

free surface 
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C. RESULTS 



A representative set of results for the first-order and 
second-order force coefficients have been generated for a 
water depth of h = 5.0, submergence depths of d = 1.5, 2.0, 
2.5, and 3.0, with a ranging from 0.15 to 1.2. All of the 
numerical results presented herein were based on 24 sub- 
divisions on the circular cylinder and 64 subdivisions per 
second-order wave length on the free surface integration. 

The surface integral was carried out from x = -L to +L since 
this was found to be adequate for convergence. Since the 
second-order wave length is half the first-order wave length, 
L, the total interval, is equal to four second-order wave 
lengths . 



The first-order horizontal and vertical force coefficients 
are presented in Figs. (3) and (4), respectively. It may be 
noted that, in general, the forces decrease with depth of 
submergence according to expectation since the wave action 
dies out with depth. 

The second-order horizontal and vertical force coefficients 
are presented in Figs. (5) and (6), respectively. Generally, 
the results show the second-order effect to be relatively more 
important as the cylinder approaches the free surface. It 
might be suspected from this that the second-order contribution 
would be much more significant in the case of surface piercing 
or floating bodies. 
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The horizontal and vertical steady-state force coefficients 
are shown in Figs. (7) and (8) , respectively. These results 
show that the horizontal steady-state force coefficient is very 
small in general. The horizontal force can be shown to be 
proportional to the momentum flux of the reflected wave and 
the reflected wave is small. In fact, in the case of infinite 
water depth. Dean [Ref. 1] showed that the reflection was 
exactly zero and accordingly the steady-state horizontal force 
is zero. 

It may be noted that the steady-state vertical force is 
positive. This results from the fact that the velocities are 
larger on the top of the cylinder and accordingly, the 
pressures reduced. 

The phase shift angles of the first-order and second-order 
forces are shown in Figs. (9) - (12) . 

To demonstrate the effect that the second-order terms have 
on the forces on the cylinder and their respective phase shift 
angles, a comparison was made between the first- and second- 
order results. E’igures (13) - (18) are plots of the horizontal 
and vertical dimensionless forces on the cylinder versus time 
over one complete wave cycle for both the first-order solution 
alone and the combined first-order and second-order effects. 
Additionally, a plot of the incident wave is included to 
demonstrate the phase shift involved (the amplitude of the 
incident wave has no significance) . A mean water depth, h, of 
5.0, a cylinder depth, d, of 1.5 and a wave height, H, of 
0.5 were used, with the wave number, a, assigned three values, 
0.25, 0.5 and 1.0. 
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Figure 3: FIRST-ORDER HORIZONTAL WAVE FORCE COEFFICIENT 
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Figure 4: FIRST-ORDER VERTICAL WAVE FORCE COEFFICIENT; 
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Figurs 5: SECOND-ORDER HORIZONTAL WAVE FORCE COEFFICIENT; 
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Figure 6: SECOND-ORDER VERTICAL WAVE FORCE COEFFICIENT 
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Figure 7: STEADY-STATE HORIZONTAL WAVE FORCE COEFFICIENT 
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Figure 8: STEADY-STATE VERTICAL WAVE FORCE COEFFICIENT 
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Figure 9: FIRST-ORDER HORIZONTAL PHASE SHIFT ANGLE 
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Figure 11; SECOND-ORDER HORIZONTAL PHASE SHIFT ANGLE; 
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Figure 12: SECOND-ORDER VERTICAL PHASE SHIFT ANGLE 
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Figure 13: HORIZONTAL WAVE FORCE; a = 0.23, h = 3.0 
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Figure 14: VERTICAL WAVE FORCE; a = 0,25 
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Figure 15: HORIZONTAL WAVE FORCE 
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Figure 16: VERTICAL WAVE FORCE; a * 0.5, h = 5.0 
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Figure 17: HORIZONTAL WAVE FORCE 
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Figure 18: VERTICAL WAVE FORCE; a = 1.0, h = 5.0 



V. CONCLUSIONS 



A computer program has been developed to carry out the 
numerical solution of the second-order wave — cylinder 
interaction problem. The first-order, second-order, and 
steady-state force coefficients were determined for the 
submerged horizontal cylinder. 

An approximate method for dealing with the second-order, 
nonhomogeneous free surface boundary condition was developed 
which appears to converge except at small values of a, i.e., 
large wave lengths. 
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// EXEC FGRTHCLG»REGI0N=350K 
//FORT . SYS I N DD * 



THIS PROGRAM CALCULATES WAVE FORCES AND THEIR PHASE 
SHIFT ANGLES FOR WAVE INTERACTION WITH A SUBMERGED 
HORIZONTAL CYLINDER 

A = 2*P I*CYL INDER RADIUS/WAVE LENGTH = SIGMA**2* 
CYLINDER RAOIUS/ACCELERATION OF GRAVITY 
H = MEAN WATER DEPTH/CYLINDER RADIUS 
D = CYLINDER DEPTH/CYLINDER RADIUS 

NPTS = NUMBER OF CYLINDER SURFACE ELEMENTS AND NODAL 
POINTS FOR NUMERICAL EVALUATION 

NSPTS = NUMBER OF FREE SURFACE ELEMENTS AND NODAL 
POINTS FOR NUMERICAL EVALUATION 



COMPLEX GI J ,GNI J ,GNJ I ,GX I J, GXJ I ,GYI J, GY J I , GI JEXT, GYY, 

1 DET f SUM l, SUM 2 , SUM 3 , SUM4 , U 1 1 S,U1ISX,U1ISY,U1ISYY, 
2U1SS,U1 SSXtUISSY »U1SSYY 
COMPLEX ALPHA (24, 24) , BET A ( 24 , 24 ) , BET AX ( 2 4 , 24) , 

1BETAY( 24, 24) ,HH( 24, 1) ,F( 24, 1) , PK(24, 1) ,F1(24,1) , 
2FS(500) ,U1 ( 24) ,U1X( 24) ,U1Y( 24) ,U2(24) ,U^SC1(24) , 
3U2SCNI ( 24) ,U2SCO ( 24) 

DIMENSION X<24) ,Y( 24) ,ANX( 24) , A NY ( 24 ) , CHY ( 2 5 ) ,CHY2( 25) 
1 , SHY ( 25) , SHY2 ( 25 ) ,COEFG(2GO) ,COEFG2( 200 ) , 

2C0SAMU ( 200,25) ,C0SAM2( 200,25) , S I NAMU ( 2 JO ,25 ) , 

3 S I NAM 2 ( 20 0,25) , AMU( 2C0 ) ,AMU4(20 0) ,SH2Yl 24) ,CH2Y( 24) 
DIMENSION Cl (2) ,C2(2) ,C3(2) , PHASE 1(2) ,PHASE2( 2) 

COMMON/ CP X /HH , PK 
C CMMON/ VAR /X , Y , A NX , ANY 

COMMCN/G SHY/C hY ,CHY2,SHY, SH Y 2 , SH2Y , CH2Y 
CCMMON/GMU/COSAMU , C0SAM2 , S I NAMU , S I NAM2 , AMU , AMU4 , COE FG, 
1C0EFG2 

C OMMON/ CON ST/H,D,DELTHE, SMI N, NPTS, NSPTS 
CCMMGN/CF I / P I 
EQUIVALENCE ( Hh( 1 ) ,F( 1 ) ) 

EQUIVALENCE ( PK ( 1 ) , F 1 ( 1 ) ) 

PI=3. 141592 



READ INPUT DATA AND CALCULATE REQUIRED REPEATING 
DATA AND ARRAYS 



CALL GEODAT ( A, A2 , ANU, ANU4, SH2AH , SH2AH2 , SHAH, SHAH2 , 
1CHAH, CH AH2 , AO, AA-, 3B , CC , DD , A02 , A A2 , BB2 , CC 2 , CD2 ) 

DO 100 1 = 1, NPTS 
DO 100 J = 1 , 1 
XV=X( I ) 

YV= Y ( I ) 

XC= X ( J) 

YC=Y( J ) 

ANX I = ANX ( I ) 

ANY I = ANY ( I ) 

ANXJ= ANX( J ) 

ANY J = ANY( J ) 

I F ( AB S ( X ( I)-X( J) ) .LT.SMIN) GO TO 50 
SHY I = SH Y ( I ) 

CHY I = CH Y ( I ) 

SHY J= SH Y( J ) 

CHY J=CHY ( J ) 

CALL GREENS ( A , ANU , S H2 AH, SH AH, CH AH, COS AMU , S INAMU , AMU , 
1C0EFG , SHY I , CHY I , SHY J , CHY J , I , J , XV , YV , XC. , YC, ANX I , ANY I , 
2ANXJ, ANYJ ,GI J ,GNI J,GNJ I , GX I J , GX J I , GY I J , G Y J I , G YY ) 

GO TO 90 
50 CONTINUE 

CALL GREEN ( A , ANU , SH2 AH , SHA H,CH AH , AO , A A , BB , CC , DD , 1 , J , 
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1XV,YV,XC, YC, ANXI , ANY I , G I J , GN I J , GX I J , G Y I J , GY Y ) 

IF(I.EQ.J) GN J I =GN I J 
IF (I.EQ.J) GX J I =GX I J 
IF (I.EQ.J) GY J I = GY I J 
I F ( I.EQ.J) GO TO 90 

CALL GREEN ( A , ANU , SH2 AH , SHA H,CHAH , AO, A A , BB , CC , DD , J , I , 
1XC,YC,XV,YV,ANXJ,ANYJ,GI JEXT,GNJI ,GXJI ,GYJI ,GYY) 

90 CONTINUE 



CALCULATE THE FIRST-ORDER ALPHA AND BETA MATRICES 



ALPHA ( I , J)=( l./PI ) *SN I J*DEL THE 
ALPHA { J , I ) = ( 1 ./PI ) *GNJI*DELTHE 
BETA (I , J)=( l./(2.*PI) )*GIJ*DELTHE 
B ETA ( J , I J = BETA ( I ,J) 

BETAX( I , J)=( l./(2.*Pl ) )*GXI J*DELTHE 
BETAX( J , I )=( I ./( 2.*PI ) USXJUDELTHE 
BETAY(I,J)=(l./(2.*Pl) )*GYI J*DELTHE 
BET AY (J, I )=(1./(2.*PI ) ) *GYJ UDELTHE 
100 CONTINUE 

DO 120 1=1 , NPTS 

ALPHA ( I , I ) =ALPH A ( I , I ) +CMPLX ( 1.0,0 .0) 

BETAX ( I , I )=BETAX ( I , I) +ANX ( I ) *C MPLX ( 0 . 5 , 0 . 0 ) 
6ETAY( I , I) =6ETAY( I , I 1 +ANY ( I )*CMPLX(0.5, 0.0) 
120 CONTINUE 



GENERATION OF THE FIRST-ORDER DISTRIBUTION FUNCTION, 
F ( I , 1 ) , BY INVERSION OF ALPHA! I , J )*F( I , 1 ) = HH( 1,1) 



CALL COMAT ( 24 , 1 , ALPHA , HH , D E T , I ND I C A i 



THE HH (1,1) MATRIX IS REPLACED BY THE DISTRIBUTION 
FUNCTION, F( 1,1) 

COMPUTATION OF THE FIRST-ORDER SCATTERING POTENTIAL 
FUNCTION, U1SC(I), AND ITS X AND Y PARTIAL DERIVATIVES 
UISCX(I) AND UISCY(I), AND COMBINATION WITH THE 
INCIDENT POTENTIAL COUNTERPARTS TO COMPUTE THE TOTAL 
POTENTIAL FUNCTION AND ITS X AND Y DERIVATIVES 



DO 155 1=1 ,NPTS 
SUM1=(0. 0,0.0) 

SUM2=( 0 .0,0.0) 

SUM3=( 0 .0,0.0) 

DO 150 J= 1 »N PT S 

SUM 1= SUNUP ( J, 1 ) *BETA( I , J ) 

SUM2=SUM2 + F( J,1 ) *BE7 AX ( I , J ) 

SUM3=SUM3+F( J , 1 UB ETAY ( I , J ) 

150 CONTINUE 

Ul( I ) =SUM1 - ( CHY ( I ) / ( A*CHAH ) )*CEXP(CMPLX(0.0,A*X( I ) ) ) 
U 1 X l I 1 =SUM2-CMPLX( 0 . 0 , 1 . 0 ) *CHY ( I ) *CEXP( CMPLX( 0.0, 

1 A*X ( I ) ) J/CHAH 

U 1 Y ( I ) = SUM3-SHY( I ) *C EX P ( CM PLX ( 0 . 0 , A* X ( I ) ) ) /CHAH 
155 CONTINUE 



EVALUATION OF THE FREE SURFACE PRESSURE DISTRIBUTION 
SOURCE STRENGTH FUNCTION, FS(L) 



DELX=0.015625*PI/A 
WRITE (6,180) DE LX 
180 FORMAT ( 5X//5X , • DELX = ' , FI 0. 5/// ) 
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SP=0.0 

CST=. 5*DELTHE/PI 
I COUNT = 1 

DU 250 L=1 i NSPTS 
SUM1 = ( 0 .0,0.0) 

SUM2=( 0 .0,0.0) 

SUM3=( 0.0,0. 0) 

SJMA-= ( 0 .0 , 0. 0) 

DO 245 1 = 1 , N PTS 

XV=SP 

YV=0.0 

XC=X( I ) 

V C=Y ( I ) 

ANX I =0 • 0 
ANYI = 1 . 0 
ANXJ=ANX( I ) 

ANY J = AN Y ( I ) 

SHYI=SHY{ 25) 

C.HYI=CHY( 25) 

SHY J = SH Y( I ) 

CHY J = CHY ( I ) 

IF (ABS (XV-XC) .LE.SMIN) GO TO 240 

CALL GREENS ( A , ANU, SH2 AH, SHAH, CHAH,COSAMU»S INAMU, AMU, 
1C0EFG , S HY I , CHY I , SHY J , CHY J , 2 5 , I , XV , YV , XC , YC , ANX I , ANY I , 
2ANXJ, ANYJ.GI J,GNI J,GNJI ,GXI J,GXJ I , GYI J , GY J I , GYY ) 

GO TO 241 

240 CONTINUE 

CALL GREEN { A , ANU , SH2 AH , SHA H , CH AH, AO , AA , BB , CC , DD , 2 5, I 
1 XV , YV , X C , YC , ANX I , ANYI ,G I J , GN1 J , GX I J , GY I J , GY Y ) 

241 CONTINUE 



SUM1=SUM1 +CST*GI J*F( 1,1) 

S i in ^ r ii m . r r r >/ r i j, r i r i \ 

UI'U 1” OUl'!£ * U O i A i J'*T i i 9 1. / 

SUM3=SUM3*CST*GY I J*F (1,1) 

SUM4=SUM4+CST *GYY*F ( 1,1) 

245 CONTINUE 

U1IS=<-1./A)*CEXP(CMPLXO.O,A*SP) ) 
U1ISX=CMPLX<0.0,-1.0)*CEXP<CMPLX(G.0,A*SP) ) 

U 1 1 SY = - TANH( A*H) *CEX?(CMPLX (0.0 , A*SP) ) 

U II SYY =— A*CE X P(CMPLX(0.0,A*SP) ) 

U1SS=SUM1 
U 1 SSX = S UM2 
U 1 SSY = S UM3 
U IS SYY= SUM4 

FS ( L ) = ( 2.*A/ 13.* ANU ) ) * ( U 1 1 S*U1 S SYY +U 1 I S Y Y*U1 SS + Ul S S’ 
1U1SSYY-6.*U1 ISY*U1SSY-3.*U1SSY*U1SSY-4.*UI ISX*U1SSX 
2-2 . *U1 S SX* U1 SSX ) 

IF ( I COUNT . E Q . 0 ) GO TO 248 
SP=FLOAT(L+l)*DELX/2. 

I C0UNT = 0 
GO TO 250 
248 SP=-SP 
I COUNT = 1 
250 CONTINUE 



, 



CALCULATION OF THE PARTICULAR SOLUTION PORTION OF THE 
SECOND-ORDER SCATTERING POTENTIAL AND ITS NORMAL 
DERIVATIVE, U2SCKI) AND U2SCNKI) 



CST=.5*DELX/PI 
DO 350 I = 1 , N PTS 
SUM 1=(0. 0,0.0) 
SUM2=(0 .0,0.0) 
SP=0. 0 
I COUNT = 1 

DO 345 L=l, NSPTS 
XV=X( I ) 
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Y V = Y( I ) 

XS = SP 

Y S=0 . 0 

ANX I = ANX ( I ) 

ANYI =ANY ( I ) 

ANX J=0 . 0 
ANY J= 1 . 0 
SHY I = SHY ( I ) 

CHYI =CHY ( I ) 

SHYJ = SHY< 25) 

CHYJ = CHY( 25) 

IF (ABSIXV-XS) .LE.SMIN) GO TO 340 

CALL GREENS ( A2 , ANU4 , SH2AH2 , SHA H2 ,CHAH2 , C0SAM2* SI NAM2, 
1 A MU4, CO E FG2 , SHY I ,CHYI , SHYJ,CHYJ, I , 25 » XV » YV . XS , YS » 

2 ANX I , ANYI, ANXJ, ANYJ,GI J ,GNI J,GNJ I »GXI J,GXJI ,GYI J, GY J I , 
3GYY ) 

GO TO 341 

340 CONTINUE 

CALL GREEN ( A2 , ANU4 , SH2 AH2 , SHAH2 , CHAH2 , A 02 , AA2 , BB2 , CC2 
1, DD2,I , 25,XV»YV,XS, VS, ANX I , ANY I , G I J , GN1 J , GX I J , GY I J , 
2GYY ) 

341 CONTINUE 
SUM1=SUM1+CST*GI J*FS( L) 

SUM2=SUM2+CST*GNI J*FS( L) 

IF ( ICOUNT.EQ.O) GO TO 290 
SP=FLOAT( L+l )*DELX/2. 

I COUNT= 0 
GO TO 345 
290 SP=-SP 

I COUNT= 1 
345 CONTINUE 

U2SC1J I ) =SUf1 1 
U2SCN1 ( I ) -SUM2 
350 CONTINUE 



CALCULATION OF THE HOMOGENEOUS SOLUTION PORTION 
OF THE SECOND-ORDER SCATTERING POTENTIAL, U2SC0I I) 



DO 500 1 = 1 , NPT S 
DO 500 J=1 , I 
XV=X< I ) 

Y V = Y( I ) 

XC = X( J ) 

YC=Y( J ) 

ANX I = ANX ( I ) 

A NY 1 = ANYI I ) 

ANXJ=ANX( J ) 

ANY J= ANY ( J ) 

IFIABSIXI I )-X( J) ) .LT.SMIN) GO TO 450 
SHY I = SH Y ( I ) 

CHY I=CHY I I ) 

SHY J=SHY ( J ) 

CHY J=CHY( J ) 

CALL GREENS I A2 , ANU4, SH2AH2 , SHAH2 ,CHAH2 , COS AM2, S I NAM2, 
1AMU4,C0EFG2, SHY I ,CHYI , SHY J , C HY J , I , J , XV , Y V , XC , YC , ANX I , 

2 ANY I , ANX J, ANY J, GI J ,GNI J,GNJ I , GX I J , GX J I , GY I J , GY J l , G YY ) 

GO TO 490 
450 CONTINUE 

CALL GREEN I A2 , ANU4 , SH2 AH2 , SHAH2 , CHAH2 , A 02 , AA 2 ,BB 2 , CC2 
1 » DD2, I , J,XV,YV, XC ,YC, ANX I , ANY I , G I J , GN I J , GX I J , GY I J , 
2GYYI 

IF(I.EQ.J) GN J I =GN I J 
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I F ( I . EQ . J ) GO TO 490 



CALL GREEN ( A2 , ANU4, S H2 AH2 , SHAH 2 , CHAH2 , A 02 , AA 2 , BB 2 , CC2 
1 » CD2 , J » 1 » XC * YC , X V » Y V » ANX J » A NY J , G l J EXT * GN J I , GXJi , GYJI , 
2GYY) 

490 CONTINUE 

ALPHA ( I,J)=!1./PI) *GY I J*DELTHE 
ALPHA! J , I ) =( 1 ./P I )* GY J 1*0 EL THE 
BETA! I , J) =< 1 ./ ( 2 .*PI ) )*GI J*DELTHE 
BETA! J, 1 ) = BE TA ( 1 ,J) 

500 CONTINUE 

DO 510 1=1 , N P T S 
PK(I,1)=2*(PK(I,1) — U2SCN1 ! I ) ) 

510 CONTINUE 

DO 520 1=1 ,NPTS 

ALPHA! I , I ) = AL PH A { I , D+CMPLX! 1.0, 0.0) 

520 CONTINUE 

CALL COMAT ( 24 , 1 , AL PHA , PK , D ET, I ND I CA ) 

DO 540 I =1 , NPTS 
SUM=! 0.0, 0.9) 

DO 530 J = 1 , N PTS 

SUM=SUM + F 1 ! J , 1 )*BETA( I , J) 

530 CONTINUE 

U2SC0 ( I ) =SUM 
540 CONTINUE 



EVALUATION OF THE TOTAL SECOND-ORDER POTENTIAL, U 2 ( I ) 



DO 550 1=1 , NPTS 

U2( I )=U2SC1! I ) +U2SC0! I ) — ( C H Y ! I ) / ! 2.* A* < SH AH**4J ) )* 
ICEXPtCMPLXIO.C, 2.*A*X! I ) ) ) 

550 CONTINUE 

WRITE (6,560) 

560 FORMAT ( 5X// 9X , • I ' , 1 5X , ' U 1 ( I ) ' , 24X , ' U1 X ( I ) ' , 2 5X , 

1 * U 1 Y ! I ) ' , 24X , ’ U 2 ! I )'//) 

00 580 1=1, NPTS 

WRITE ( 6,570) I , U 1 ( I ) , U IX ( I ) , U 1 Y ( I ) , U2 < I ) 

570 FORMAT ( 5X , I 5 , 4 ( F 1 6 . 6 , F 14 . 6 ) ) 

580 CONTINUE 



EVALUATION OF THE FIRST-ORDER, SECOND-ORDER PERIODIC, 
AND SECOND-ORDER STEADY STATE FORCE COEFFICIENTS AND 
THE PERIODIC FORCE PHASE SHIFT ANGLES 



SUM1=(0. 0,0.0) 

SUM2=(0 .0,0.0) 

SUM3=( 0 .0,0.0) 

SUM4= ( 0 . 0 , 0. 0 ) 

SUM5=0 « 0 

SUM6=0 . 0 

DO 720 1=1, NPTS 

SUM1 = S UM 1 +U 1 ( I)*ANXm#DFLTHE 

SUM2= SUMZ + U1 ! I) *ANY! I ) *DELTHE 

SUM3 = SUM3 + ( ( 6.*ANU*ANU*U2 ( I ) /A ) -U1X ! I ) *U1 X ( IJ-U1Y! I) 
1*U1Y! I ) )*ANX! I ) *DELTHE 

SUM4 = SUM4 + ( ( 6. * ANU* ANU*U2 ( I ) /A ) -UiX! I ) *U IX ( I) -U1Y ( I ) 
1*U1Y( I ) )*ANY< I ) ^D EL T HE 

SUM5=SUM5+ ( CABS !U IX! I ) *U1 X ! I ) M-C ABS ( U 1 Y ! I )*U1Y( I ) ) 

1 + ANU*ANU/ A** 2- 1 .0)*ANX( I ) *DELTHE 
SUM6=SUM6+- ( CABS ( U1X ( I ) *U1 X ( I ) ) +CABS { U1Y ( I ) *U1 Y< I ) ) 

1 + ANU*ANU/A**2-1 .0)*ANY( i )*DELTHE 
720 CONTINUE 

C1{1)=CABS(SUM1*A) 

C 1 (2) =CABS ( SUM2*A ) 

C 2 ( 1 ) = CABS ! SUM3* A* A/ ( 4 . *ANU ) ) 
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C2(2)=CABS(SUM4*A*A/(4.*ANU) ) 

C3( 1) =SUM5*A*A/< 4. *ANU) 

C3( 2 )=SUM6$A*A/( 4 . *ANU) 

PHASE 1 ( 1 ) = A T AN2 ( A I M AG ( S JM 1 * A J ,REAL(SUML*A)) 

PHASE 1 ( 2 ) = A T AN2 ( A I MA G ( SUM2 * A ) , RE AL ( SUM2* A J ) 

PHASE 2 ( 1 ) = AT AN2 ( A I MAG ( SUM3*A*A/ ( 4 . *ANU ) ) , RE AL ( SUM 3* A*A 
1/ (4.*ANU) ) ) 

PHASE2 ( 2)=AT AN2( AIMAGt SUM4*A*A/ (4.*ANU) ) , RE AL( SUM4*A*A 
1/(4. *ANU) ) ) 

WRITE (6,730) C 1 ( 1) , C 1 ( 2 ) , C 2 ( 1 ) , 02 ( 2 ) , C 3 ( 1 ) ,C 3 ( 2 ) , 

1 PHASE 1 ( 1 ) , PHASE1 ( 2 ) , PHASE 2 ( 1 ) , PHASE2( 2) 

730 FORMAT (5X///5X, 'CL< 1 ) = ' , F8 . 5 , 5X , * C 1 ( 2) = ' , F 8. 5/ / 5X , 
l'C2(l)=',F8.5,5X,'C2(2)= , ,F8.5//5X, , C3(l)=',F8.5,5X, 

2' C3(2)=* ,F8.5//5X, 'PHASEl ( 1)=' ,F8.5,5X, • PHASE1(2)=' , 
3F8.5//5X, ' PHA SE 2 ( 1)=',F8.5,5X,'PHASE2(2)=',F8.5//) 

STOP 

END 



THIS SUBROUTINE READS THE INPUT GEOMETRICAL DATA AND 
CALCULATES THE FI RSI 200 ROOTS OF AMU* I AN ( AMU*H) + 

ANU = 0 AND AMU4*TAN(AMU4*H) f ANU4 =0. IT ALSO 
GENERATES ARRAY'S OF CERTAIN FUNCTIONS AND COEFFICIENTS 
USED IN GREEN AND GREENS AS WELL AS THE HH AND PK 
MATRICES 



SUBROUTINE GEODAT ( A , A2 , ANU , ANU4 , SH2 AH , SH2AH2 , SHA H , 

1 SHAH2 , CMAh ,CHAH2 , AO, A A, BB , CC , DD , A02 , AA2 , B6 2 , CC2 , DD2 ) 

COMPLEX HH ( 24, 1 ) , PK ( 24 , 1 ) 

DIMENSION X ( 24 ) , Y ( 24 ) , ANX ( 24 ) , ANY ( 24 ) ,CHY( 25) ,CHY2(25) 
1,SHY(25), SHY 2(25) , CC EFG ( 200 ) , C0EFG2 ( 200 ) , AMU ( 200 ) , 
2AMU4( 200) ,C0SAMJ ( 20C , 25) , COSAM2 ( 200 , 25) , S I NAMU( 200 ,25 ) 
3 , SINAM2 (200,25) ,SH2Y( 24) , CH2Y( 24) ,XX( 24) 

C OMMON / C PX / HH , P < 

CCMMON/VAR/X, Y, ANX , ANY 

COMMON /C’SHY/’CHY , CHY2 , SHY , SHY2 , SH2Y , CH2Y 
C OMMON / C- MU/ C OS A M U , C 0 S A M2 , S I NA MU , S I N AM2 , A MU , AMU4 , C OE F G , 
1C-0EFG2 

COMMON/C ON ST /H, D , DEL.THE , SMI N,NPTS , NS PTS 
COMMON/C.PI/PI 

RE AO (5, 5) A , H, D, SM IN, NPTS , NSPTS 
5 FORMAT (4F 10, 5, 214) 

ANU=A*TANH< A*H) 

SH2AH=S INH<2.*A*H) 

CHAH=COSH( A*H) 

SHAH=SINH( A*H) 

AO=T AN H( A*H) 

AA=1./CHAH**2 
E B=-H*SHAH/( CHAH**2) 

CC=-H*H*( 1 .-SHAH**2)/ (3.*CHAH**3) 

DD = H**3=i SHAH*( 2 . *CHAH**2 + 3. *( I . - SHAH** 2 ) ) / ( CH AH**4* 1 2 ) 

DELTHE=2.*PI/NPTS 

THETA=DELTHE/2 . 

DO 15 1=1, NPTS 
X ( I ) =COS ( THETA ) 

Y ( I ) = S I N( THE T A ) -D 
ANX ( I ) = X ( I ) 

ANY ( I ) = Y ( I )+D 

SHY( I) = SINH( A*(H*Y(I) ) ) 

C HY ( I ) = COSh( A*( H+Y ( I ) ) ) 

SH2Y ( I) =SINH(2.*A*(K*Y( I ) ) ) 

C H2Y ( I)=C0SH(2.*A*<H*-Y(I)J) 

HH( It 1 ) = ( 2 ,/CHAH)*CMPLX( ANY ( I ) *S HY ( I ) , A NX ( I )*CHY< I ) )* 
1CEXP( CMPLX(0.0 , A*X(I))) 

PK( 1 , 1 ) =( ANY ( I)*Sri2Y(I )+CMPLX(0 .0, l. 0) *ANX ( I ) *CH2Y( I ) ) 
1*CEXP(CMPLX(0.0,2.0*A*X(1 )) )/SHAH**4 
THETA=THETA+DELTHE 
15 CONTINUE 

SHY (25) = SH AH 
CHY ( 25 ) =ChAH 
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B=ANU*H 
DO 50 K= 1 » 200 
XX( 1)=PI*K 
CO 25 1=1,20 
I I = I 

Y Y=XX( I ) 

XX ( 1 + 1) = XX( 1 )-ATAN2(3 ,YY) 

IF(ABS( < X X ( I + 1 ) - XX ( I i ) /XX( I +1) ) .LT..0001 ) GO TO 30 
25 CONTINUE 

IF (II. GE. 20 1 GO TO 27 
GO TO 30 

27 WRITE (6,28) 

28 FORMAT (5X,'AMU ROOT DOES NOT CONVERGE'//) 

30 CONTINUE 

AMU(K) =XX (II ) /H 

COEFG (K) = 2.*PI*( AMU(X) *AMU( K) +ANU*ANU) / ( ANU*AMU( K) -H 
1 AMU ( K ) * ( AMU ( K ) * AMU ( K ) +ANU-ANU ) ) 

NPTS1=NPTS+1 

DO 40 I=1,NPTS1 

IF (I .EC.NPTS1 ) GO TO 35 

COS AMU ( K, I )=COS( AMU(K)*(H+Y( I ) ) ) 

S I NAMU ( K , I )=SIN( AMU(<)=MH + Y( I ) ) ) 

GO TO 4C 

35 C OS AMU ( K , I ) =COS ( AMU ( K ) *H) 

S I NAMU ( K, I )=S IN ( AMU(K)*H) 

40 CONTINUE 
50 CONTINUE 

ANU4=4.*AMU 
B2=ANU4*H 
DO 70 K = 1 , 200 
XXII )=PI*K 
DO 55 1=1, 20 



.T..G001 ) GO TO 60 



1 1 = I 

Y Y=XX ( I ) 

XX(I+1)=XX(1)-ATAN2(82,YY) 

I F ( AS S ( ( XX ( I i 1)~XX( 

55 CONTINUE 

IF ( II .GE.20) GO TO 57 
GO TO 6C 

57 WRITE (6,58) 

58 FORMAT (5X,'AMU4 ROOT DOES NOT CONVERGE*//) 

60 CONTINUE 

AMU4(K) =XX ( I I )/H 

C0EFG2 ( K ) =2. *P I * ( AMU4 ( K) *AMU4( K) +ANU4*ANU4) /( ANU4* 

I AMU4 ( K ) -H*AMU4 ( K ) * ( AMU4 ( K ) *AMU4 ( K ) +ANU4* ANU 4 ) ) 
NPTS1=NFTS+1 
DO 65 I = I , NPTS 1 
IF (I.EQ.NPTS1) GO TO 68 
C0SAM2 ( K, I )=COS( uMU4(K)*(H+Y( I ) ) ) 

SINAM2 (K, I )=SI N( AMU4 ( K )* ( H + Y( I ) ) ) 

GO TO 65 

68 C0SAM2 ( K , I ) =COS ( AMU4 ( K ) *H) 

S I NAM2 ( K, I) = SIN(AMU4(K)*H) 

65 CONTINUE 
70 CONTINUE 
XX( 1 )=ANU4 
DO 80 1=1,20 
I I = I 

Y 2 = XX ( I ) 

XX ( 1+1 ) =4. *ANU/TANH(Y2*H) 

IF ( A3S ( t XX( 1+1 ) -XX( I )) /XX( 1+1 )) .LT. 0.0001) GO TO 85 
80 CONTINUE 

WRITE (6,32) 

62 FORMAT (5X,'A2 ROOT DOES NOT CONVERGE'//) 

85 CONTINUE 
A 2=XX ( I I ) 

SH2AH2=SINH( 2.*A2*H) 

CHAH2=C0SH( A2*H) 

SHAH2= S INH( A2*H) 

A02=TANH( A2*H) 

AA2=1. /CHAh2**2 

Bt3 2=-H*SHAH2/( C. HAH 2**2 ) 
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CC2=-H*H*( l.-SHAH2**2)/(3. *CHAH2**3 ) 

DD2=(H**3)*SHAH2* ( 2 . *CHAH 2** 2 + 3 . *< 1 . -S H AH2**2 ) ) / 

1 ( (CHAH2**4)*12. ) 

DO 90 1=1 ,NPTS 
S HY 2 ( I ) =S IMH( A2* ( ri+Y( 1))) 

CHY2( I )=C0SH(A2*(H+Y( I ) ) ) 

90 CONTINUE 

WRITE (6,95) A, A2 ,D,H ,ANU,ANU4, NPTS, SMIN,NSPTS 
95 FORMAT ( 5X / / / 3 X , • A = • , F 1 0 . 5 , 4X , • A 2 = • , F 1 0 . 5 , 4 X , * D= » , 
1F10.5»4X» ' H= • , F 1 0 » 5 , 4 X , ' ANU = ' , F10.5, 4X, • ANU4=* ,F10. 5/ 
24X, • NPTS= 5 , 13, 4X , 1 SMIN=* , F10.5,4X, , NSPTS = ' , 13/// ) 
RETURN 
END 



THIS SUBROUTINE CALCULATES G AND ITS DERIVATIVES, 

GX , GY , GN , AND 3 Y Y , 3Y USE OF THE INTEGRAL FORM FOR THE 
CASE (X(I) - X(J)) LESS THAN SMIN 



SUBROUTINE GREEN ( A , A NJ , S H2 AH, SH AH ,CHAH , AO , AA , BB , CC , DD 
1,I,J,X,Y,XI , ETA, ANX,ANY,G,GN,GX,GY,GYY) 



COMPLEX G, GN 
COMPLEX GX ? GY, GY Y 

DIMENSION TEST! 200) , TESTT( 100) , TESTTT( 100) ,SUMOX( 15) , 
1 SUMOY ( 15) ,NNN< 15) , TEST YY( 200) 

C CMMON/ CON ST/H,D,DELTHE, SMIN, NPTS, NS PTS 
COMMON /CP I /P I 

P1(Y,ETA,XX,AMU)=C0SH(AMU*(Y+H) )*COSH(AMU*< ETA+H) )* 

1 COS ( AMU*XX )/ (COSH ( AMJ*H)**2 ) 



P2X(Y,ETA,X,XI , AMU) =-AMU*COSH( A MU* (ETA+H ) ) *SIN( AMU* 

1 ( X-X I ) )*COSH( AMU* (Y+H) ) / ( COSH ( AMU*H ) ** 2 ) 
P2Y(Y,ETA,X,XI,AM'J)=AMU*C0SH( AMU* ( ETA + H ) ) *S INH( AMU* 

1 (Y+H) ) *COS ( AMU* ( X-X I ) ) / ( C OS H( AMU*H ) **2 > 

P2YY ( Y , C T A , X , X I , A MU ) = A MU* AMU *C OSH ( AM J* ( Y +H ) ) *CCS H ( AMU* 
1 ( ETA + H) )*COS( AMU*( X-aI ) ) / ( COSH < AMU*H ) **2 ) 



Q1 ( Y, ETA, XX, AMJ) =-Pl ( Y , ET A, XX , AMU) * ( AMU-A)/ (AMU* 

IT ANH ( AMU*H ) - AN J ) 

Q2X(Y , ETA, X, XI, A MU )=-P2X(Y, ET A , X , X I , AMU ) * ( A MU- A) /( AMU* 
1 T ANH( AMU*H)-ANJ ) 

Q2Y(Y,ETA,X,XI»4MU)=-P2Y(Y»ETA»X»XI,AMU)*( AMU- A) /(AMU* 
IT ANH( AMU*H) -ANU ) 

Q2 YY ( Y , F TA , X , X I , AMU) =-P2YY ( Y , ET A , X , X I , AMU ) * ( AMU- A ) / 

1 ( AMU*T ANH( AMU*H ) -ANU) 

Q10( Y ,ETA,XX)=-COSH( A*( Y + H) ) *COS H ( A* ( ET A + H ) )*COS ( A*XX) 
1* A/ ( COS H( A *H ) ** 2 * ( ( A-A-ANU*ANU) *H + ANU) ) 

Q20X ( Y , ETA,X ,XI ) =A*C3SH( A*( ETA+H) ) *S IN( A* (X-X I ) ) *COSH 
1 ( A*( Y + H ) ) *A/COSH( A*H)**2*( ( A* A- ANU* AN J ) *H+ ANU) ) 
Q20Y(Y,ETA,X,XI ) =-A*COSH ( A* ( ETA+H ) ) * S I N H ( A* ( Y +H ) ) *COS 
1 ( A* ( X-X 1 ) ) *A/ ( COSH ( A*H ) **2*( l A* A- ANU* ANU ) *H+ANU) ) 

Q20YY ( Y , ETA, X,XI ) =-A*A*COSH ( A* ( Y+H) ) *COS H( A*( ETA+H) ) * 
1C0S(A*( X-XI ) )*A/(C0SH(A*H)**2*( ( A* A- ANU* ANU ) *H+ANU ) ) 
Q1S< Y , ETA, XX ,AMU ) =-Pl( Y , ETA, XX, AMU) / ( A0+ AMU*H* ( AA +BB* 

1 ( AMU- A) +CC*( AMU-A)**2+DD*( A MU- A ) **3) ) 

Q2XS( Y , E T A , X , X I , A.-iU ) =-P2X ( Y , ETA, X ,XI , AMU )/( AO + AMU*H* 

1 ( AA+BB* ( AMU- A) + C C " ( AMU- A ) **2+DD*( AMU-A) **3) ) 

Q2YS( Y, ETA ,X, XI , AMU ) =-P 2 Y ( Y , E T A , X , X I , AMU) / ( AO+AMU*H* 

1 ( AA+BB* ( AMU- A ) +CC*( AMU- A ) **2 + DD* ( AMU- A) **3 ) ) 



Q2YYS ( Y , E T A, X , X I ,AMU)=-P2YY(Y,ETA,X,XI, AMU ) / ( A0+ AMU*H* 
1 ( AA + BB* ( AMU-A ) +CC*( AMU-A ) **2+DD* ( AMU- A ) **3) ) 

FUN1 (Y , ET A, XX, A MU) =EX P ( -A MU *H ) *S I NH ( AMU* ET A)*SINH 
1 ( AMU*Y ) *COS ( AMU* XX ) / ( AMU* COSEK AMU*H ) ) 

FUN2( Y , ETA, XX) =4 .*3. 14159*C0SH( A* (Y+H) ) *COSH( A* 

1( ETA + H) )*COS( A*XX ) /( 2 . * A * H + S H 2 A H ) 

FUNR ( X , Y , X 1 , ET A , a.NX , ANY ) = ( ( X-X I ) *ANX + ( Y + ET A ) *ANY ) / 
l( (X-XI)**2+(Y+ETA)**2) 

FUN3X(Y , E T A, X , X I , AMU ) =EXP ( -AMU*H ) *S I NH < AMU* ETA) *S IMH 
1 ( AMU*Y ) * S I N ( AMU* (X-XI) ) /C OS H ( AMU *H ) 

FUN BY { Y , ETA, X, XI , AMU) =-EXP ( -AMU*H ) * S I NH ( AMU*ETA ) *COSH 
1 ( AMU*Y ) *COS( AMU* ( X-X l ) ) /C OS H ( A MU*H ) 

F UN3YY ( Y , E TA , X , X I , AMU ) =-E X P ( - AMU*H) *S I NH ( AMU* ET A ) * 
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IS INH( AMU*Y)*AMU*AMU*COSI AMU*(X-XI ) )/COSH ( AMU*H) 
FUN4I Y,ETA,XX,ANX,ANY,SIGN)=4.*3. 14 1 59* A*COSH ( A* 

1 ( ETA + H ) ) *< ANY*SINHI A*( Y+H ) ) *C3 S I A*XX ) -ANX*COSH(A# 
2(Y+H))*SIN( A*XX)*SIGN)/(2.*A*H+SH2AH) 



EVALUATION OF THE FINITE INTEGRAL IN G TO DETERMINE 
THE SIZE OF THE SUBDIVISIONS 



XX = ABS ( X— X I ) 

IFIX.LT.XI) SIGN =-1.0 
IF(X.GT.XI) S I GN = 1 . 0 
IF(X.EO.XI) S I GN=0 . 0 
DO 50 N = 1 , 15 
DELMU=2*A/ I o*N+3 ) 

AMU=0. 0 
SUM=0 . 0 

F0=(Q1< Y, ETA, XX, AMU ) -Q 1 0 ( Y , E T A , XX) )/( AMU- A) 

LL = 6*N + 3 
DO 40 NN=1 » LL 

I F( ABS( AMU-A ) .LT . .00001 ) GO TO 10 
IF (A3S( AMU+DELMU-A ) .LT. .00001) GO TO 10 
1FIABSI AMU+DELMU/3.-A) .LT. .00001 ) GO TO 10 
IFIABSl AMU+2.*QE LMU/3 .-A) .LT. .00001) GJ TO 10 
F 1= < Q 1 ( Y, ETA, XX, AMU+DELMU/3. I -QIC! Y, ETA, XX) )/( AMUf 
1DELMU/3.-A) 

F 2= ( Q1 ( Y f ETA,XX,AMU*2.*DELMU/3. )-Q10(Y,ETA,XX))/ 

1( AMU+2.*DELMU/3.-A) 

F3= I Q1 I Y, ETA, XX, AMU+OE LiMU )-Q10(Y , ET A » XX ) )/( AMU+DELMU 
1-AI 

GO TO 30 

10 F 0= ( Q1 ( Y, ETA, XX, AMU+DELMU) - Q10 < Y , ETA, XX) ) /{ AMUfDELMU 
1 — A ) 

GO TO 40 

SUM= IDELMU/3. ) * I FO +3 .*F H 3. *F2 *F3 ) {-SUM 

F 0=F 3 

AMU=AMU+DELMU 
TEST(N)=SUM 
IF(N-l) 50,50,45 
MN=N-1 
MM=6*MN+3 

1F(ABS( (TEST(N)-TEST(N-l) )/TEST(N)).LT..010) GO TO 60 
CONTINUE 
CONTINUE 
PV1F=2.*SUM 



EVALUATION OF THE FINITE INTEGRAL IN GN USING 
2* A/MM SUBDIVISION SIZE 



D£LMU=2*A/MM 
AMU=0. 0 
SUMX=0 . 0 
S UMY = 0 . 0 

F 0= ( Q2X ( Y , ET A , X , X I , AMU ) -Q2 OX ( Y , ETA, X, XI ) )/( AMU-A) 
Y0=(Q2Y (Y»ETA»X,XI , AMU ) -020Y I Y , ET A , X , X I ) )/( AMU-A) 

DO 80 NN= 1 , MM 

IFIABSl AMU-A) .LT. .00001 ) GO TO 70 

IFIABSl AMU+DELMU-A) .LT. .00001) GO TO 70 

IF I ABS ( AMU+-DELMU/3 .-A ) .LT. . 00001 ) GO TO 70 

IFIABSl AMU+2.<OELMU/3. -A) .LT. .00001) GO TO 70 

F1 = (Q2X (Y, ETA, X, XI , AMU+-DELMU/3 . ) -Q2QX ( Y , ETA,X,XI ) ) / 

1 I AMU* DE LMU/3 .-A ) 

F2=(Q2X(Y,ETA,X,XI , AMU*2 . *D L LMU / 3 . ) -Q20X I Y , ET A , X , X I ) )/ 
1 I AMU+-2. *DELMU/3 .-A) 

F3=(Q2X(Y ,ETA,X»XI , AMU *DE LMU ) -Q2 OX I Y , ET A , X , X I ) )/ 

1 I AMU+-DE LMU-A ) 

Y1=(Q2Y (Y,ETA,X,XI , AMU+DE LMU/3 . )-Q20Y(Y , ETA,X,XI) )/ 

1 l AMU+DELMU/3.-A ) 

Y2= I Q2Y I Y , ET A , X , X 1 , AMU + 2 .*DE LMU/3 . ) -Q20Y I Y , ET A, X, X I ) )/ 



30 

40 

45 

50 

60 
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1( AMU+2. *DELMU/3.-A) 

Y 3= ( Q2Y ( Y i ETA.X, X l , AMU +DELMU ) ^Q20Y ( Y , ET A , X , X I ))/ 
1 ( AMU+DELMU-A ) 

GO TO 75 
70 CONTINUE 

E 0= ( Q2X ( Y , ETA, X, XI , AMU + OEL MU ) - 02CX ( Y , E T A , X , X I ) )/ 
1 ( AMU+DELMU-A ) 

Y0= ( Q2Y (Y, ETA, X,XI , AMU+DELMU ) -Q20Y (Y , ET A , X , XI ) )/ 
1 ( AMU+DELMU-A) 

GO TO 30 
75 CONTINUE 

SUMX= SUMX + ( DELMU/8. ) * ( F 0 +3 . * F 1+3 . *F2 + F 3 ) 

SUMY= SUMY + (DELMU/8. )* ( YO +3 . * Y 1 + 3 . *Y2 +Y3 ) 

F 0 = F3 
Y0 = Y3 

80 AMU= AMU+DEL MU 
PV2FX=2 .*SUMX 
PV2FY=2 .*SUMY 



EVALUATION CF THE FINITE INTEGRAL IN GY Y USING 
2* A/MM SUBDIVISION SIZE 



DELMU=2*A/MM 
AMU=0 . 0 
SUMYY = 0.0 



YY0=(Q2YY(Y,ETA,X,XI , AMU) -Q20Y Y ( Y , ETA , X , X I ) )/(AMU-A) 

DO 350 NN = 1 5 MM 

IF(ABS( AMU-A) .LT. 0.00001) GOTO 320 
IF (ABS ( AMU + DELMU-A ) .LT. 0.00001 ) GO TO 320 
IF ( ASS ( AMU + DELMU/3. -A) .LT. 0.00001 ) GO TO 320 
IF (ABS(AMU+2.*DcLMU/3.-A) .LT. 0.00001) GO TO 320 
YY1 = ( Q2YY ( Y, ETA,X,XI , AMU + DELMU/3 * J-Q20YY ( Y, ETA, X, XI))/ 
1 (AMU+DELMU/3. -A) 

YY2= ( Q2YY ( Y , ETA ,X , XI , AMU + 2 . *DE L MU/3 . i -320YY (Y,ETA, X, 

IX I ) )/( AMU + 2.*QELMU/3.-A) 

YY3=( Q2YY( Y, ETA, X,XI , AMU + DELMU ) -Q20YY ( Y , ETA,X,XI ) )/ 

1 ( AMU+DELMU-A) 

GO TO 325 
320 CONTINUE 

YY0=(Q2YY( Y, ETA, X, XI , AMU + DELMU ) -Q20YY ( Y , ETA,X,XI ) )/ 

1 ( AMU+DELMU-A ) 

GO TO 350 
325 CONTINUE 

SUMYY = SUMY Y+ (DELMU/8. ) * ( Y YO+3 . * YY 1 +3 . * Y Y 2+Y Y3 ) 

Y YO=YY 3 

350 AMU=AMU + DE LMU 
PV2FYY=2.*SUMYY 



EVALUATION OF THE INFINITE INTEGRAL IN G, GN , AND GYY 
SIMULTANEOUSLY 



AMU=2*A 
DELMUO= DELMU 

F0=Q1(Y ,£TA,XX,AMU)/( AMU-A) 

FOX=32X(Y, ETA.X, XI , AMU) /(AMU-A) 

F0Y=Q2Y(Y,ETA,X,XI ,AMU) /( AMU-A) 

F0YY=Q2YY(Y,ETA,X,XI ,AMU) /( AMU-A) 

SUM=0 . 0 
SUMX=0 . 0 
SUMY=0. 0 
SUMYY=0.0 
DO 100 NN= 1 ,200 
DO 95 N = 1 , 20 

F1=Q1 (Y, ETA, XX, AMU+DELMU/3. )/( AMU+DELM'J/ 3 .-A) 

F 1X=Q2X( Y, ETA ,X, XI , AMU+DELMU/3 . ) / ( AMU + DELMU/ A .-A > 

F 1 Y=Q 2 Y ( Y , ETA , X , X I , AMU + DELMU/3 . i / ( AMU + DELMU/3. -A ) 
F1YY = Q2YY (Y,LTA,X,Xl, AMU + DELMU/3. ) / ( AMJ + DE L MU/3 . - A ) 
F2=Q1 ( Y , ETA, XX, AMU+2. *DELMU/3 . ) / ( AMU + 2 . * DEL MU/ 3 . - A ) 
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F2X=Q2X (Y , ETA,X, XI , AMU + 2 . *DELMU/ 3 . )/(AMU+2.#DELMU/3. 
1-A) 

F2Y=Q2Y(Y, ETA,X, XI , AMU* 2 . *DE LMU/3 . ) / < AMU + 2 . *DELMU/ 3 . 

I - A ) 

F 2YY=Q2YY ( Y , ETA, X, XI , AMU* 2 . *DE LMU/3 . I / ( AMU*- 2. *OEL MU/ 3. 
1 - A ) 

F3 = Q1 ( Y , ETA, XX ,AMU«- DE LMU ) / ( AMU+DELM J - A ) 

F 3X=Q2X ( Y , ETA ,X ,XI ,AMU + OELMU) / ( AMU+DELMU-A) 
F3Y=U2Y(Y,ETA,X» X I , AM U+DELMU)/ < AMU+DELMU- A ) 
F3YY=Q2YY(Y,ETA,X»XI, AMU+DE LMU ) / < AMU + DE L MU- A) 

SUM= (DELMU/3. ) * ( FO +3 . *F H- 3. * F 2 +F3 ) + S DM 
S UMX = <DELMU/8.)*<F0X+3.*F1X+3.*F2X+F3XH-SUMX 

SUMY= (DcLMU/3.)*(FOY + 3.*FlY4-3.*F2Y + F3YH-SUMY 
SUMYY=$UMYY+ OELMU/8. ) * < FO Y Y + 3 . * F 1 Y Y + 3 . * F2 Y Y + F3 Y Y ) 

FO = F3 
FOX=F3X 
FOY=F3Y 
F OYY=F3 YY 

A SM= EXP ( AMU* ( E T A + Y ) )/AMU 
I F ( ASM. LT . .3001 I GO TO 86 
95 AMU= AMU + OE LMU 
86 TEST(NN)=SUM 
T ESTT ( N N ) = SUMX 
TESTTTI NN) =SUMY 
T ESTYY ( NN i =SUMY Y 
IF(NN-l) 100,100,97 

97 CONTINUE 

IFIASM.LT. .00010) GO TO 105 
IF(ABSISUM) .LT..0001) GO TO 98 

I F < ASS < (TEST(NN)-TEST(NN-l) )/TEST(NN)) .GT..OOI) 

1 G 0 TO 100 

98 IF(ABS( SUMX) .LT. .0001) GO TO 99 

I F ( AB S ( ( TESTTI NN)-TESTT(NN-i) ) / TE ST T ( NN i) . GT. . 00 I ) 

.1 GO TO 100 

99 IFIABSI SUMY) .LT. .0001) GO TO 93 

I r / A o r / ' "F CTT T ' mr- c T T "F * ' * * • 1 ' ' / Tr" F t T t t i » r* -r \ 

nAuol \ jloI \ i i 1 1 w ” i c. j i i i \ i s i n — JL 4 t / ICO I I I i I \\ \ t / »UI » * UU1 / 

1 GO TO 100 

93 IF(ABS( SUMYYJ.LT. 0.0001) GO TO 96 

IF (ABS< { TESTY Y(NN) -TESTY Y INN- 1) ) /TESTY Y ( NN ) ) .GT..001) 
1 GO TO 100 
96 CONTINUE 
GO TO 105 
100 CONTINUE 
102 WRITE! 6, 103) 

103 FORMAT (3X3 4H INFINITE INTEGRAL DID NGT CONVERGE) 

105 CONTINUE 

P VI 1=2 . *SUM 
PV2 I X = 2 . *S UMX 
PV2IY = 2. * SUMY 
P V2I YY = 2 . * SUMY Y 
D ELMU= . 2/H 

IF (J.EQ.25) GO TO 240 

AMU=0 . 0 

SUM=0. 0 

SUMX=0 . 0 

SUMY=0 . 0 

SUMYY=0 . 0 

F 0=0.0 

FOX=FUN3X(Y,ETA,X,XI , AMU) 

F0Y=FUN3Y(Y,ETA,X,XI ,AMU) 

F0YY=FUN3YY( Y, ETA ,X,XI ,AMU) 



DO 200 NN = 1 ,100 

IF( AMU*H.GT.5. ) DELMU=.l 

DO 195 N= 1 , 20 

F IX = FUN3X( Y ,ETA,X,XI , AMU*DE LMU/3 . ) 
F2X=FUN3X( Y ,ETA,X,XI , AMU*- 2 . *DE LMU/3 . ) 
F3X=FUN3X ( Y , ETA, X ,X I , AMU* DE LMU ) 
F1Y=FUN3Y(Y,ETA,X,XI , AMU+DE LMU/3 . ) 
F2Y=FUN3Y< Y,ETA,X,XI , AMU* 2 . *DE LMU/3. ) 
F3Y=FUN3Y( Y,ETA,X,XI , AFiU+DE LMU ) 
F1YY=FUN3YY(Y,ETA,X,XI , AMU+DcLMU/3 . ) 
F2YY=FUN3YY(Y»ETA,X,XI , AM‘J + 2 . *DE LMU/3 . ) 
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F3YY=FUN3YY( Y,ETA,X,XI , AMU+DE LMU) 

I F( AMU+DE LHU.LT • .001) GO TO 120 
FI- FUN 1 ( Y , ETA, XX, AMU + DE LMU/ 3. ) 

F 2=FUN 1 (Y,cTA,XX, AMU+2.*DELMU/3. ) 

F3=FUN 1 (Y , ETA, XX , AMU+DE LMU ) 

GO TO 130 
120 CONTINUE 

F 1 = EXP ( - ( AMU + DEL MU/ 3 . )*H) *C0S( ( AMU+OELMU / 3 . )*XX) * 

1 ( AMU+DE LMU/ 3. ) * Y *E T A/ COSH ( { AMU + DE LMU/ 3. ) *H) 

F2=EXP ( -( AMU+2. *DELMU/3 . ) *H ) -COS ( ( AMU +2 . *DELMU/3 . )*XX> 
1*< AMU+2.-DELMU/3 . ) *Y*ETA/C0SH( ( AMU+2 . *CE LMU /3 . ) *H ) 
F3=EXP(-( AMU+DE LMU) *H)*CO$( ( AMU * DELMU ) *X X )*( AMU + DE LMU) 
1*Y*ETA/CGSH< ( AMU+DE LMU ) *H ) 

130 CONTINUE 

SUM= (DELMU/3. ) * ( FO +3 .*F 1 +3 . *F2 +F3 ) +SUM 

SUMX* ( DELMU/8 .)* ( FO X+3 . * F IX +3 . *FZX +F 3 X ) + SUMX 

SUMY= (DELMU/8. ) * ( FOY+3 . *F 1 Y + 3 . *F 2Y+F 3 Y ) + SUMY 

SUMYY=SUMYY+ ( DELMU/8. ) * ( F OY Y+ 3 . *F 1 YY+ 3 . *F2Y Y+F3YY ) 

F0=F3 

F0X=F3X 

F 0 Y = F 3 Y 

F0YY=F3 YY 

ASM=EXP( AMU*( ETA+Y) ) 

IF( ASM .LT. .0001 J GO TO 196 

195 AMU=AMU +DELMU 

196 TE ST ( NN ) = SUM 
TESTT( N N ) = SUMX 
TES TTT ( NN ) =SUMY 
TESTYY ( NN) =SUMYY 
IF(NN-l) 200,200,199 

199 CONTINUE 

IF (ASM. LT. .00010) GO TO 205 
I F(ABS( SUM) .LT. .0001 ) GO TO 206 

I F ( AB S ( ( TEST ( NN) -TEST (NN-1) ) / TES T ( NN ) ) . GT . . 00 10 ) 

1 GO TO 200 

rs a I c / Aor / muv \ it o ^ ^ T ' ^ r> to 00-7 
C JO il (MUJl UUilA J *L l c iJuui I L 7 U t CUI 

IF (Ad SI ( TESTT(NN)-TESTT (NN-i) )/TESTT (NN) ) .GT..00I0) 

1 GO TO 200 

207 IF ( ABS ( SUMY) .LT. .0001 ) GO TO 204 

I F ( ABS ( ( TESTTT( NN) -TES TTT (NN-1) ) / TESTTT ( NN ) ) . GT . . 001 0) 
1 GO TO 200 

204 IF(A8S( SUMYYJ.LT. 0.0001) GO TO 203 

IF ( ABS ( ( T ESTYY ( NN )- TESTYY ( NN- 1 ) )/( TESTYY ( NN) ) ) . 

1GT. .001 ) GO TO 200 

208 CONTINUE 
GO TO 205 

200 CONTINUE 
WRITE ( 6 , 202) 

202 FORMAT ( 3X14HN0 CONVERGENCE) 

205 CONTINUE 
GINF=-2*SUM 
GXINF=2 .*SUMX 
GY I NF = 2 . *SUMY 
GYYINF=2. *SUMYY 
GNS I NG= .5 

IF (I.EQ.25) GO TO 218 
I F{ l.EQ . J) GO TO 220 
A I J J- I - J 

THETA1=A'3S< AIJJ) *DELTHE 

AINJJ=I +NPTS-J 

T HET A2 = A8S ( A IN J J ) *DELTHE 

AJNII=I-J-NPTS 

T HE T A3 = ABS ( A JNI I ) *DELTHE 

THETA=THcT A1 

IF(THETA2.LT. THETA) THETA = THET A2 
IF ( THETA3.LT. THETA) THETA = THETA 3 
IF(THETA.GT..15) GO TO 218 

GS I NG=D fc LT HE *ALOG( 2. ) + 2 . * ( - DEL T HE / 2 . + ( D EL TH E/ 4 . + 

1 THETA/ 2 . )*AL0G( T HET A/ 2 . +DELTHE/4 . ) - ( THE T A / 2 .-DELTHE / 4 . 
2 ) >--4003 ( THETA/2. -DELTHE/4. )-( DELTHE/4. + TH ETA/2. )** 3/18. 
3+ ( THETA/2. -DEL THE /4. ) **3/18 .-( THE TA/2. +D ELTHE/4 . ) **5/ 
4< 180. *5. ) +( THE T A/ 2. -DEL THE/ 4. )**5/< 160. *5. ) -(THETA/2. + 
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5DELTHE/4.)**7/( 2835. *7 . ) + ( THETA/2. -DELTHE/4. I **7 / 
6(2835. *7. ) ) 

GXSI NG = £ X. — X I )/( < X-XI ) ** 2 +- ( Y-ET A ) **2 ) 

GYSING=( Y-ETA)/< (X-XI ) **2 + ( Y-E TA ) **2) 

GO TO 230 

218 GSING=DELTHE*ALOG( SQRT( ( X-X I ) **2 * ( Y-E T A ) ** 2 ) ) 

GXSING= ( X-XI )/< (X-XI)**2t( Y-ETA)**2) 

GYSING= (Y-ET A) /( ( X-XI ) **2 + ( Y-E T A ) **2 ) 

GO TO 2^0 
220 CONT INUE 

GSING=DELTHE*ALOG (DEL THE/ 2. ) -D ELTH E- I DE L TH E**3 ) / 

1{ l8--5=l£>.)-(DELTHE^^5)/( 180. *5. *256. I -( 0ELTHE**7 ) / 
2(2835. *7. *256.^16. ) 

GXSING=GNSINb*ANX 
GYSING=GNSING*ANY 
230 CONTINUE 

G = GSI NG/DcLTHE- ALOG ( SQRT ( (X-XI ) **2 + ( Y + ET A )**2 ) ) +PV1F 
1«-PV1H-GINF+CMPLX(0. » - 1 . )*FUN2( Y , ETA, XX) 

GX=GXSI NG-FUNR (X,Y,XI,ETA,1.0,0.0)+PV2FXi-PV2IX + GXINF + 
1FUN4< Y , ETA f XX, 1 . 0 ,0.3, SI3N)*CMPLX (0.0,- l .0) 

GY = GYSING-FUNR( X,Y,XI,ETA,0.0,1.0) +PV2FY+PV2IY+GYINF+ 
1FUN4( Y, ETA, XX, 0.0, 1.3, SIGN) *CMPLX (0.0, -1.0) 

IF (I.LE.NPTS) 30 TO 235 

FURYY= 1 ./SORT ( ( { X-XI ) **2fr( Y-ET A ) **2 ) **3 I - ( ( Y-ET A 1**2 + 
1( Y-E TAJ )/SQRT( ( (X-XI J **2M Y-ET A ) **2)**5 ) 

FUNYY= 1 . /SQRT( ( ( X-XI )**2M Y + ETA )**2) **3 ) -( < Y+ETA)**2+ 

1 < Y+ETA) ) /SORT ( ( (X-XI ) **2M Y +ET A) **2) **5 ) 

GYY=FURYY+ FUNYY+ PV2FYY+P V2 I YY* GYY INF +CMPLX ( 0.0,- 1 .0)* 
1FUN2< Y, ETA,XX)*A*A 
235 CONTINUE 

IF (I.EQ.25) GO TO 250 

GN=GNS 1 NG-FUNR ( X , Y , X I , ETA, ANX, ANY )+{ PV2FX+PV2 IX+GXINF) 
1* ANX+ ( PV2FY*PV2IY+GYINF)*ANY+FUN4(Y,ETA,XX, ANX, ANY, 
2SIGN)*CMPLX(0.0,-1.0) 

GO TC 250 
240 CONTINUE 

G=- ( P V 1 F + PV 1 I )+CMPLX(0. ,-l. )*FUN2(Y,ETA,XX) 

GN=- ( PV2FX+PV2IX) *ANX- ( P V2F Y + P V2 I Y ) *ANY * FUN4( Y,ETA, XX, 
IAN X, ANY, SIGN) * C M P L X ( 0 . 0 , -1.0) 

250 CONTINUE 
RETURN 
END 



THIS SUBROUTINE CALCULATES G AND ITS DERIVATIVES, 

GX, GY, GN , AND GYY , BY USE OF THF SERIES FORM FOR THE 
CASE (X(I) - X(JJ) GREATER THAN SMIN 



SUBROUTINE GREENS ( A , ANU, SH2 AH , SHAH ,CHAH , COSA MU , S I N AMU 
1 , AMU, COEFG , SHY I , CHY I , SHY J ,CHY J , I , J,XV, Y V ,XC ,YC, ANXI , 

2 ANY I , ANXJ , ANYJ , GI J ,GN I J ,GNJ I , »XI J ,GX J I , GYI J ,GY JI , GYY ) 

COMPLEX GI J , GN I J , GN J I , GX I J , GX J I , C-Y I J , GY J I , GYY 
DIMENSION C 0 S A M U ( 200,25) , S I N AM U ( 200, 25) ,ANU(200) , 
1C0EFG ( 2 OC ) 

DIMENSION TEST1(4G) ,TEST2(40J ,TEST3(40) ,TEST30(40) , 

1 TEST4( <*0> 

C GMMON/ CONST /H, D,DELTHE, SMI N,NPTS, NS PTS 

CCMMON/C PI /PI 

SUM 1=0 . 0 

SUM2=0. 0 

SUM3=0 . 0 

SUM30=0 .0 

S UM4=0 . 0 

DO 20 N = 1 , 40 

DO 15 K K= 1 ,5 

K={ N-l ) *5 + KK 

SN I =S INAMU ( K , I ) 

SNJ=SINAMU (K, J ) 

C S I = COS AMU ( K , I ) 

CSJ=COSAMU(K, J) 
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AK= AMU ( K) 

CCK=COEFG ( K ) 

VAL=AK*ABS(XV-XC ) 

IF ( VAL .GE. 75.0) GO TO 10 
EXI J=EXP(-VAL ) 

GO TO 1 1 

10 E X I J = 0 . 0 

11 CONTINUE 

SUM1 = SUM1 +COK S I *C S J^EX I J 
SUM2=SUM2+CCK*CSI *CSJ*EXI J*AK 
SUM3=SUM3*C0K*AK*SNI*CSJ*EX I J 
SUM30=SUM30*C0K*AK*SNJ*CS1*EXI J 
SUM4=SUM4+ AK*AK*COK*C SI -CSJ-EXI J 

15 CONTINUE 

TE ST 1 I N ) = SUM 1 
TEST2I N )=SUM2 
TEST3 ( N ) =SUM3 
TEST30IN) =SUM30 
T EST4 1 N ) = SUM4 
IFIN-l) 20,20,16 

16 IFIABSI SUM1) .LE.O.OOOOOl) GO TO 17 

IFIABSI ( TEST1(NJ-TEST1(N-1) ) / TEST 1( N ) J.GT..0010) GO 
1 TO 20 

17 IFIABSI SUM30) .LE .0.000001) GO TO 18 

IFIABSI I TEST30IN )-TEST30I N-l ) ) /TEST3 DIN) ) .GT. .0010) 

IGO TO 20 

18 IFIABSI SUM3) .LE. 0.000001) GO TO 19 

IFIABSI I TEST3IN)-TEST3IN-1) J/TEST3IN) ) .GT..0010) GO 
1T0 20 

19 IF (ABSI SUM2) .LE .0.000001 ) GO TO 21 

IFIABSI < TEST2IN)- TEST 2IN-1) J/TEST2 IN) ) .GT.. 00 10) GO 
1T0 20 

21 IFIABSI SUM4) .LE.O.OOOOOl) GOTO 30 

IFIABSI (TEST4(N)-TEST4(N-1))/TE$T4(N)).GT.. 0010) GO 
1 TO 20 
GO TO 30 

20 CONTINUE 

WR I TE I 6 , 2 5 ) I,J 

25 FORMAT! 10X23HGREENS DID NOT CONVERGE , 2X2HI = 12 , 

1 2X2HJ= 12) 

30 CONTINUE 

IF IXV.LT.XC) S I GN=- 1.0 
IF (XV.GT.XC) SIGN =1.0 
IF IXV.EQ.XC) SIGN =0.0 

TERM4=4.*PI*CHYI *CHYJ*SIN( A*ABS( XV-XC ) ) / 1 2. *A*Hf SH2AH) 
T ER M 5= 4 . * P I * C H Y I *C H Y J *C 0 S I A*ABS! XV-XC) ) / I 2 . *A*H + SH2 AH ) 
TERM6=A*TERM5*SIGN 
TERM7=A*TERM4*SI GN 

TERMS=4.*P I*A*SHYI*CHYJ*COS( A* ABSI XV-XC ) ) /( 2.*A*H* 
1SH2AH) 

TERM9=4.*PI*A*SHYI*CHYJ*SIN I A* ABSI XV-XC ) ) /I 2.*A*H + 
1SH2AH) 

TERM80=4.*PI*A*SHYJ*CHYI*C0SI A* ABSI XV-XC ) ) / I 2 .*A*H + 
1SH2AH) 

TERM90=4.*PI*A*SHY J*CHYI*SIN( A*ABS(XV-XC ) )/ (2.*A*H + 
1SH2AH) 

TERM10=A*A*TERM4 
TERM11=A*A*TERM5 
IF IJ.EQ.25) GO TO 70 
G I J=C MP LX I SUM1 + TE RM4 , -T ER M5 ) 

GX I J= CMPLXI-SIGN*SUM2<-TERM6,TERM7) 

GY I J = CMPLX (-SUM3+-TERM9 , -TERMS) 

IF I I . E 0 . 2 5 ) GO TO 40 

GX JI=CMPLX( SIGN*SUM2-TERM6 ,-TERM7) 

GYJI=CMPLX(-SUM30+TERM90, -TERM80) 

GNI J = CMPLX( -ANX I*S IGN--SUM2 + A NX i *TERM6-A NY I * SUM3* ANY I * 

1 T E R M 9 , A N X I # T E R M 7 - A N Y I * T E k M 8 ) 

GNJI=CMPLX( ANXJ*SI GN* S UM2 - ANX J * T ERM6- AN Y J*SUM30+ANY J* 
1TERM90, -ANXJ*TERM7-ANYJ*TERM80) 

GO TO 60 
40 CONTINUE 

GYY=CMPLXI -SUM4+TERM10 ,-TERMl 1 ) 
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60 CONTINUE 
GO TO 30 

70 GI J=CMPIX<-SUMI-TE RM4 , -TERM 5) 

GNI J=CMPLX (ANXI-S I GN* SUM2- A NX I * T ERM6 + AN Y I*SUM3-ANYI* 
1TERM9, ANXI*TERM7-ANYI*TERM8 ) 

80 CONTINUE 
RETURN 
END 



THIS SUBROUTINE INVERTS COMPLEX MATRICES TO SOLVE 
THE MATRIX EQUATION AL PHA< I , J ) *F ( I , 1 ) = HH (1,1) AND 
ALPHA(I f J)*F1( I i 1) = PK(I,li 



SUBROUTINE COMA T ( N , M , A , 6 , D , I ) 

INTEGER C,H,R,Q, Z 

COMPLEX A,8,D,TT,P 

D I ME NS I ON A ( N, N ) , 6 ( N, M ) , C ( 1 00 , 3 ) 

D = ( 1 .0,0.0) 

DO 20 J = i , N 

20 C ( J , 3 ) = 0 

DO 2L K = 1,N 
TT = (0.0,0.0) 

T « 0.0 
DO 4 J = 1 , N 

IF { C t J , 3 ) .EQ. 1) GO TO 4 

DO 5 H = I ,N 

IF ( C ( H , 3 ) - 1) 15,5,12 

15 IF (T .GE. CABS(A(J,H) )) GO TO 5 

R = J 

Q = H 

T = CABS! A ( J , H) ) 

5 CONTINUE 
4 CONTINUE 

C < Q , 3 ) = C I Q, 3 ) + i 
C ( K , 1 ) = R 

C ( K , 2 ) = Q 

IF (R . EQ. Q) GO TO 11 
D = -D 

DO 6 L = 1 ,N 
TT = A ( R , L ) 

A I R , L ) = A I Q , L ) 

8 A ( Q, L) = TT 

I F ( M .LE. 0) GO TO 11 
DO 2 L = 1 ,M 
TT = B ( R , L ) 

B(R,L) = B ( Q , L ) 

2 B(Q,L) = TT- 
11 P = A(Q,Q) 

A(Q,Q) = (1.0, 0.0) 

DO 13 L = 1 , N 
13 A ( Q , L ) = A ( 0 , L ) / P 

IF (M .LE. 0) GO TO 1 
DO 3 L = 1 , M 

3 B(Q,L) = B(Q,L)/P 
1 DO 21 Z = 1 , N 

IF ( Z . EQ. Q) GO TO 21 
TT = A ( Z , Q ) 

A(Z,Q)=( 0.0, 0.0) 

DO 16 L = 1,N 

16 A ( Z , L ) = A ( Z , L ) - A ( Q , L ) *TT 
IF (M o LE. 0) GO TO 21 

DO 17 L = 1 , M 

17 B(Z,L) = B(Z,L) - B ( Q, L ) ^TT 

21 CONTINUE 

DO 19 II = 1 , N 
L = N * 1 - I I 

I F ( C ( L , 1 ) .EQ. C ( L , 2 ) ) GO TO 19 
R = C ( L , 1 ) 

Q = C( L , 2 ) 
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DO 7 K = 1 » N 
TT = A ( K » R ) 

A(KtR) = A ( K t Q) 

7 A l K t Q ) = TT 
19 CONTINUE 

DO 18 K = 1 * N 

IF ( C ( K » 3 ) .NE. 1) GO TO 12 
18 CONTINUE 
I = 1 
50 RETURN 
12 I = 2 
GO TO 50 
END 



//GO.SYSIN DD * 
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